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Abstract 


The  discrete  ordinates  method  is  widely  used  to  solve  the  Boltzmann 
transport  equation  for  neutral  particle  transport  for  many  engineering 
applications.  Source  iteration  is  used  to  solve  the  discrete  ordinates  system  of 
equations,  but  can  be  slow  to  converge  in  highly  scattering  problems.  Synthetic 
acceleration  techniques  have  been  developed  to  address  this  shortcoming; 
however,  recent  research  has  shown  synthetic  acceleration  to  lose  effectiveness  or 
diverge  for  certain  problems. 

LTC  Wager  introduced  an  alternative  to  source  iteration  and 
demonstrated  it  in  slab  geometry.  Here  the  method  is  further  developed, 
enhancing  efficiency  in  various  ways,  and  demonstrated  in  XY-geometry  as  well 
as  slab  geometry.  It  is  shown  to  be  efficient  even  for  those  problems  for  which 
diffusion-synthetic  and  transport-synthetic  accelerations  fail  or  are  ineffective. 
The  method  has  significant  advantages  for  massively-parallel  implementations. 
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Distribution  Iteration:  A  Robust  Alternative  to  Source  Iteration 
for  Solving  the  Discrete  Ordinates  Radiation  Transport  Equations 

In  Slab  and  XY  -  Geometries 

I.  Introduction 

The  time-independent,  single  energy  group,  linearized  Boltzmann 
Transport.  Equation  (BTE)  for  non-multiplying  systems  can  be  written: 

[n.v+<7((r)Mr,n)  =  Jrfa/a-J(r,^.n)Kr,no+^(i:,a)s  (1.1) 

where  y/  is  the  angular  flux;  <JI  is  the  total  cross  section;  (X  is  the  scattering 
cross  section;  and  qext is  the  external  source  (5:  2).  The  BTE  in  this  form  is  an 
integro-differential  equation  that  is  coupled  in  space  and  angle.  The  discrete 
ordinates  method  discretizes  the  BTE  in  space  and  angle  and  the  resulting 
system  of  equations  is  widely  used  for  solutions  to  the  BTE  for  many  engineering 
applications. 

This  research  demonstrates  a  new  method  that  is  a  robust,  flexible  and 
rapid  way  of  solving  the  discrete  ordinates  system  of  equations.  Various 
techniques  have  been  applied  to  solve  the  discrete  ordinates  equations  with 
varying  degrees  of  success.  A  brief  review  of  several  techniques  follows. 
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A.  Background 


1.  Source  Iteration 

One  technique  that  is  commonly  used  to  solve  the  discrete  ordinates  system 
of  equations  is  known  as  source  iteration  (SI).  The  BTE  can  be  written  in  an 
operator  notation: 

L^  =  S  y/  +  E,  (1.2) 

where  L^  =  [Q-  V +  <T,(r)fy/r(f ,Q)  is  the  streaming  and  collision  operator; 

=  Jr/Q'cr  (r,  Q'  ■  Cl)i//(r,  Q!  ’)  is  the  scattering  operator  or  (within  group) 
scattering  source,  and  E  is  the  emission  source1  (which  includes  scatter  into  the 
group  from  other  groups  in  a  multigroup  formulation).  As  an  iterative  scheme, 

SI  is  written  (5:  2): 

lV/+1)=S^(/)+A.  (1.3) 

The  BTE  is  discretized  in  angle  and  space.  An  initial  estimate  of  the 
scattering  source  is  made.  The  right  side  of  equation  (1.3)  is  treated  as  the 
source  for  this  method,  the  sum  of  both  the  scattered  particles,  as  determined  by 
the  integral  and  the  emission  sources  in  the  material.  The  discretization  in  angle 
allows  the  integral  for  the  scattering  source  to  be  evaluated  using  a  quadrature 
rule  with  the  initial  flux  estimate  for  N  directions  to  determine  the  source.  If  the 
initial  guess  for  the  scattering  source  is  0,  then  the  l- th  iteration,  or  estimate  of 

1  My  notation,  S  and  E  ,  rather  than  q scat  and  qext ;  was  chosen  to  reserve  subscripts  and  superscripts 
for  other  uses. 
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the  angular  flux,  is  due  to  particles  that  have  scattered  at  most  l-l  times  (5:  2). 
The  number  of  scatters  that  must  be  modeled  determines  the  speed  with  which 
SI  converges. 

Source  iteration  has  been  used  for  many  years,  but  has  several 
shortcomings.  For  problems  that  are  dominated  by  scattering  with  little  or  no 
absorption,  the  SI  method  may  take  many  iterations  to  converge  and  require 
impractical  compute  times.  Further,  in  highly  scattering  problems,  the  difference 
between  two  iterations  may  meet  the  convergence  tolerance  before  the  true 
solution  is  reached.  This  phenomenon  is  known  as  false  convergence.  Techniques 
to  speed  the  convergence  have  been  studied  with  varying  degrees  of  success  (18: 
1,5:  1).  Currently  there  is  no  technique  to  speed  the  convergence  that  works  for 
all  problems  in  two  dimensions. 

2.  Synthetic  Acceleration 

Methods  have  been  developed  over  the  past  40  years  to  accelerate  or 
rapidly  converge  SI,  particularly  for  diffusive  type  problems.  One  technique  that 
is  commonly  used  is  synthetic  acceleration,  which  is  at  least  a  two  stage  iteration 
scheme.  The  first  stage  is  a  normal  iteration  from  SI,  with  a  change  of  the 
iteration  subscript  from  equation  (1.3): 

Li/sa+i)  =Si/sU)  +E .  (1.4) 
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The  intent  of  the  second  stage  is  to  find  a  low  order  approximation  to  add  to 
!//  2  as  a  better  approximation  to  the  exact  solution  y/ .  Subtracting  equation 
(1.4)  from  equation  (1.2)  and  solving  for  the  exact  solution: 

yr  =  y/0^  +  (L-S)^1  S (y/u+^ -y/U))  .  (1.5) 

Finding  (L  —  S)~l  to  a  high  order  is  as  difficult  as  solving  the  original  problem; 
therefore,  a  low  order  approximation  is  used  where  M  =  (L—  S)-1  is  easier  to 
compute.  The  synthetic  acceleration  scheme  then  is: 

y/(l)  =y/<l+i)  +  MS(y/iI+{>  -y/(,))  .  (1.6) 

Diffusion  synthetic  acceleration  (DSA)  and  transport  synthetic 
acceleration  (TSA)  are  two  commonly  used  synthetic  acceleration  methods.  The 
DSA  scheme  uses  a  diffusion  approximation  as  the  low  order  approximation, 
while  the  TSA  scheme  uses  a  simplified  transport  operator,  for  example  a  smaller 
angular  quadrature,  as  the  low  order  approximation.  For  homogenous  material 
problems,  these  techniques  have  been  highly  effective  (5:  2-3). 

Adams  and  Larsen  presented  a  comprehensive  review  of  these  methods,  as 
well  as  others,  along  with  their  strengths  and  limitations  (2:  139).  For  problems 
with  severe  spatial  heterogeneities,  DSA  in  multiple  dimensions  has  been  shown 
to  degrade  significantly  and  TSA  has  been  shown  to  diverge.  Additionally,  a  new 
consistent  differencing  derivation  is  needed  for  each  new  type  of  problem  with 
DSA,  and  TSA  still  has  difficulties  for  problems  that  are  highly  scattering.  As 
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the  authors  state,  there  is  strong  interest  in  new  methods  that  are  efficient  and 


easy  to  implement. 

3.  Angle  Iteration 

Wager  developed  a  new  method  to  solve  the  BTE  that  could  be  a  practical 
replacement  for  source  iteration.  His  method  is  called  Angle  Iteration  (AI)  and 
uses  iteration  on  the  cell  edge  flux  distribution  to  rapidly  converge  on  a  flux 
solution.  His  method  does  not  converge  falsely.  His  work  showed  promising 
results  but  was  only  demonstrated  in  slab  geometry  for  isotropic  scatter. 

His  method  begins  by  treating  the  discretization  in  angle  and  space  as  a 
system  of  equations,  representing  the  flux  for  all  directions  as  a  vector  and  the 
spatial  relations  as  a  matrix  multiplying  the  flux  vector.  For  a  single  cell  (cell  i) 
in  one  dimension  for  any  spatial  method,  the  outgoing  flux,  the  incoming  flux, 
average  flux  and  average  source  relations  in  his  notation  are: 

Wont,  ~  ^Ol}fr'ml  "*"-^-05,^4  >  (l-^) 

Ya,  =  +  K aeEa  ,  (1.8) 

SA=TSiYAr  (1.9) 

In  these  equations,  KOI ,  KOSA  ,  KOEA  ,  KAI ,  KASA  ,  and  KAEA  ,  are 
diagonal  matrices  of  transport  coefficients.  Each  element  is  the  quantity  of  flux 
(out  or  average  in  the  cell)  constituted  by  the  uncollided  (first-flight)  streaming 
of  a  unit  quantity  of  flux,  scattering  source  or  emission  source.  Only  the  first 
flight  flux  is  included  in  each  K  ;  the  flux  of  scattered  particles  is  included  as  the 
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first  flight  of  the  (previously)  scattered  source  particles  (S).  The  quantities 
SA  and  Ea  are  the  variables  that  represent  the  average  scatter  and  average 
emissions  in  a  cell.  Also,  Y,s  represents  the  scattering  cross  sections  with  the 
appropriate  quadrature  weights  to  calculate  the  scattering  source  from  the  cell 
average  angular  flux  (16:  2-26). 

Equations  (1.9)  and  (1.8)  can  be  substituted  into  equation  (1.7)  to  solve  for 
the  vector  of  cell  face  fluxes  out  of  a  cell  in  terms  of  the  vector  of  incoming  fluxes 
and  the  vector  of  emission  in  the  cell: 


=  (Kw,  +K„_  ISi(I-K,S;  + 

(KOI +K»S.  XPI -K,s  Zp-'K.Jfi,, 


(1.10) 


k  OE ,  1  OS ,  ‘-‘S, 

1  : 


The  factor  (/  — Xs.)  1  in  the  above  equation  is  the  sum  of  an  infinite 
geometric  series  as  long  as  <  1 .  Further,  each  term  in  the  sum  models 

a  scattering  event  within  a  cell.  The  factor  (/  — K^.  )_1  therefore  models  all 

numbers  of  scatters  that  a  particle  can  have  before  leaving  the  cell  (16:  2-55).  I 


call  this  infinite  within-cell  scatters.  This  is  different  than  SI  which  models  each 
scattering  event  separately  and  hence  has  difficulties  with  dominantly  scattering 
problems.  Equation  (1.10)  can  be  given  in  a  compact  notation  which  represents 


the  matrices  in  the  outer  parentheses  as  a  single  matrix: 


You,,  =mOIiYini+mOEiEAl  •  (1-11) 

While  this  relation  does  model  infinite  within  cell  scatters  (15:  2-31),  it  accounts 
for  contributions  to  the  flux  from  other  cells  in  the  slab  only  indirectly,  through 
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y/in  .  Scattering  among  cells  is  addressed  by  representing  equation  (1.11)  as  a 
coupled  system  of  equations  across  all  the  cells: 

ocEt.  (1.12) 

Further,  the  incoming  flux  in  a  cell  is  the  outgoing  flux  from  adjacent  cells, 
(except  at  exterior  boundaries)  hence: 

*„  =  ?' F„.  (1-13) 

where  P  is  the  permutation  matrix  that  reorders  the  outgoing  flux  vector 
appropriately.  Substituting  equation  (1.13)  into  equation  (1.12)  yields  (after  some 
algebra): 

F.,=(P,.(I-MwP)r,M(,££).  (1.14) 

where  Pm  is  a  permutation  matrix  that  reorders  the  matrix  to  be  of  minimum 
bandwidth.  This  system  of  equations  fully  couples  angle  and  space  to  get  a  flux 
solution,  but  is  impractical  to  solve  for  fine  angular  and  spatial  resolution 
because  it  is  the  full  set  of  simultaneous  discrete  ordinates  equations  (16:  2-44). 
Wager’s  AI  method  makes  use  of  the  strengths  of  both  equations  (1.11)  and  (1.14). 

In  the  AI  method,  the  outgoing,  within  cell  flux  is  solved  using  equation 
(1.11).  The  flux  solution  is  then  collapsed  into  two  directions.  The  collapsing  is 
done  by  summing  (integrating)  the  fluxes  in  a  given  hemisphere  over  the 
hemisphere.  The  collapsed  flux  is  used  to  solve  equation  (1.14)  for  two  directions 
across  the  spatial  grid  for  what  he  called  the  global  problem.  The  new  flux 
solution  from  the  global  problem  is  apportioned  back  into  the  original  cell 
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representation  using  the  original  flux  distribution  (16:  2-73).  This  initial 
distribution  is  normalized  to  create  flux  weights.  A  flux  weight  for  a  direction  is 
the  flux  moving  in  that  direction  divided  by  the  sum  of  all  the  flux  moving  along 
the  same  hemisphere.  A  similar  flux  weight  can  be  defined  for  the  opposite 
direction,  as  well  as  the  average  flux  in  a  cell,  and  the  scattered  and  emission 
source  in  a  cell.  This  process  of  collapsing,  solving  and  apportioning  gives  a 
better  estimate  for  the  cell  edge  flux  and  can  be  used  to  solve  for  an  updated  cell 
edge  flux  (15:  2-66-69).  The  updated  cell  edge  flux  can  be  collapsed  with  new 
flux  weights  and  the  process  repeated.  This  process  describes  one  iteration.  The 
iterative  process  is  continued  until  a  convergence  tolerance  is  met  (16:  2-73). 

The  AI  method  has  been  tested  using  two  positive  spatial  methods:  step 
characteristic  (SC)  and  exponential  characteristic  (EC).  In  both  cases,  it  was 
shown  to  be  reliable  and  to  rapidly  converge  across  a  broad  range  of  cross 
sections  and  a  full  range  of  scattering  ratios  for  these  positive  spatial  methods 
(16:  6-1). 

Despite  the  success  of  the  AI  method,  there  are  several  issues  to  address: 
angular  quadrature  choices,  cell  particle  flow  variable  representation,  and 
coupling  of  the  scattering  among  cells.  These  issues  will  be  developed  and 
addressed  in  the  next  two  chapters. 

Additionally,  the  AI  method  was  demonstrated  for  spatial  quadratures 
that  only  required  the  calculation  of  a  zeroth  spatial  moment  of  the  flux  in  a  cell. 
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The  higher  spatial  moments  for  the  nonlinear  EC  method  were  found  through  a 


root  solving  routine.  Implementation  of  linear  first  spatial  moment  methods  was 
not  yet  derived.  In  addition,  the  effect  of  using  non-positive  linear  spatial 
methods  in  the  AI  method  needed  to  be  examined. 

The  AI  method  was  demonstrated  in  slab  geometry.  An  extension  to 
multiple  dimensions  required  addressing  two  issues:  how  to  incorporate  the  flux 
scattering  from  the  orthogonal  directions;  and  how  to  efficiently  communicate  cell 
information  about  cell  emissions  and  absorptions  across  the  spatial  mesh.  For 
one  dimension,  the  global  problem  resulted  in  a  pent.a-diagonal  matrix  which  can 
be  solved  efficiently.  A  similar  coupled  global  problem  in  two  dimensions  needed 
to  include  the  scattering  terms  as  well. 

B.  Motivation 

Despite  the  challenges  that  needed  to  be  addressed  for  the  AI  method,  the 
results  demonstrated  in  slab  geometry  showed  promise  that  a  flexible,  robust- 
method  could  be  developed  and  demonstrated  in  multiple  dimensions.  Further, 
Wager’s  tests  in  slab  geometry  suggested  that  this  new  method  could  overcome 
difficulties  that  SI  and  synthetic  acceleration  methods  have  for  particular 
problems  in  XY-geometry. 
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C.  Goal  of  the  Research 


The  goal  of  this  research  was  to  develop  and  demonstrate  a  new  algorithm 
for  rapid  solutions  of  the  discrete  ordinates  equations  in  two  dimensions.  It  is 
desirable  that  the  algorithm  be: 

Robust  -  able  to  handle  a  broad  range  of  cross  sections  and  scattering 
ratios  without  significant  changes  in  convergence  rates; 

Flexible  -  able  to  easily  implement  additional  spatial  methods  without 
requiring  another  derivation  and  change  to  the  algorithm.  The  method  should 
also  be  able  to  change  angular  quadratures  with  no  changes  to  the  algorithm; 

Parallelizable  -  although  the  method  was  implemented  and  demonstrated 
on  a  desktop  machine,  it  is  desirable  that  the  method  be  parallelizable  to  be  able 
to  handle  large  problems  efficiently;  and 

Readily  extendable  to  3D  -  the  methodology  used  in  deriving  and 
implementing  the  method  should  provide  a  clear  path  to  implementing  the 
algorithm  in  three  dimensions. 

D.  Objectives: 

1.  Extend  the  method  to  2-d  Cartesian  Geometry. 

2.  Use  other  spatial  and  angular  quadratures  to  inherit  correct  diffusion 

limits. 
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3.  Evaluate  the  utility  of  a  partial  current  problem  (a  finite-volume 
particle  conservation  formulation)  as  an  alternative  to  Wager’s  use  of  partial 
range  angular  integrals  of  the  directional  flux  for  coupling  cells  in  a  global 
problem. 

4.  Formulate  the  method  to  minimize  the  size  of  the  global  problem  when 
applied  to  higher  order  linear  methods. 

5.  Demonstrate  success  where  both  DSA  and  TSA  fail  or  become 
ineffective  and  extend  testing  to  even  more  challenging  problems. 

6.  Evaluate  the  ability  of  a  PARDISO-based  direct  solver  routine  (6:  11- 
1)  to  solve  the  partial  current  problem  efficiently. 

7.  a)  Maximize  the  opportunity  for  parallelization, 
b)  Enhance  serial  performance. 

8.  Distribution  iteration  should  have  the  desirable  properties  described  as 
goals  of  the  research. 

E.  Scope 

The  scope  of  this  research  is  to  derive  and  implement  a  new  method  for 
solutions  to  the  discrete  ordinates  equations  using  linear  spatial  methods  for  slab 
and  XY  -  geometry  with  discrete  ordinates  quadratures.  Slab  geometry  testing, 
for  both  zeroth  and  first  spatial  moment  methods,  was  used  to  validate  method 
choices  for  XY-geometry  testing.  Implementation  of  the  DI  method,  for  both 
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zeroth  and  first  spatial  moment  methods,  show  general  performance  of  the 
method  for  a  variety  of  parameters.  In  addition,  tests  in  XY-geomet-ry  show  the 
improvement  the  new  method  has  over  other  methods  currently  used  to  solve 
these  same  equations.  The  code  implementation  was  written  to  be  able  to 
demonstrate  this;  it  is  not  intended  to  be  incorporated  in  a  production  code. 

F.  Assumptions  and  Limitations 

This  research  uses  linear  spatial  methods  that  provide  solutions  to  the 
time  independent,  mono-energetic  BTE  for  isotropic  scatter  and  non- multiplying 
systems  in  two  dimensions.  Energy  dependence  is  not  tested  explicitly. 
Nevertheless,  the  emission  source  can  include  scatter  into  a  group  from  other 
groups,  so  the  derivations  would  apply  to  a  multigroup  formulation  without  loss 
of  generality.  Similarly,  my  testing  assumes  isotropic  scattering,  but  this 
influences  only  the  numerical  values  of  the  elements  of  the  scattering  matrix,  Z  • 
Extension  to  anisotropic  scatter  requires  only  the  formulation  of  Z  consistent 
with  the  anisotropic  scatter  approximations  to  be  employed. 

The  new  method  solves  the  discrete  ordinates  equations  and  therefore 
inherits  the  strengths  and  weaknesses  of  the  angular  and  spatial  quadratures  and 
the  cross  section  approximations  used. 
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G.  Approach 


The  first  step  was  to  examine  the  appropriate  choice  for  cell  particle  flow 
for  implementation  in  the  distribution  iteration  method.  A  change  in  the 
representation  for  the  problem  by  transforming  the  angular  flux  representation  of 
the  AI  method  into  a  current  representation  for  the  cell  transport  coefficients  is 
appropriate.  This  allows  changing  the  “global”  problem  for  the  AI  method  into  a 
partial  current  problem.  This  was  done  for  several  reasons.  This  is  a  more 
physical  problem  which  is  based  on  the  conservation  of  particles  as  opposed  to  a 
pseudo  scalar  flux  which  was  used  in  the  AI  method.  Using  this  representation, 
the  extensions  to  three  spatial  dimensions  are  more  apparent  and  the  same 
methodology  can  be  used.  Also,  test  problems  in  chapter  three  showed  that  the 
method  converges  in  fewer  iterations  for  this  representation.  The  flux  weights 
used  in  the  AI  method  are  replaced  by  current  distributions  on  the  cell  edges. 

This  motivates  the  name  of  the  new  method:  distribution  iteration  (DI). 

The  angular  integrals  described  in  the  BTE  were  done  using  an  angular 
quadrature.  The  discrete  elements  quadrature  used  in  the  AI  method  did  not 
meet  the  diffusion  limit,  which  is  needed  for  the  highly  scattering  problems  this 
research  attempted.  The  discrete  ordinates  quadratures  that  are  commonly  used 
do  meet  the  diffusion  limit.  Two  different  quadratures  were  implemented  for  two 
reasons,  to  compare  with  previously  published  results  and  to  demonstrate  the 
flexibility  of  the  method  in  implementing  different  angular  quadratures. 


13 


Different  methods  to  couple  the  scattering  among  cells  are  presented  in 


chapter  two  and  tested  for  efficiency  in  chapter  three.  Next,  the  DI  method  for 
zeroth  spatial  moment  methods  is  reviewed  and  first  spatial  moment  methods  in 
slab  geometry  are  derived  in  chapter  three.  Implementation  of  two  spatial 
methods,  step  characteristic  (SC)  and  linear  discontinuous  (LD)  are  covered. 

Test  problems  were  used  to  validate  choices  for  the  distribution  iteration  method 
implementation  in  XY  -  geometry. 

The  DI  method  is  derived  for  zeroth  and  first  spatial  moment  methods  in 
two  dimensions  in  chapter  four.  The  following  methods  were  implemented:  step 
characteristic  (SC);  weighted  diamond  difference  (WDD);  linear  characteristic 
(LC);  and  linear  discontinuous  (LD).  The  zeroth  spatial  moment  methods  (SC 
and  WDD)  validate  the  extension  from  one  dimension  to  two  dimensions,  while 
the  derivation  and  implementation  of  the  more  complicated  first  moment 
methods  (LC  and  LD)  further  demonstrate  the  flexibility  of  the  method.  The 
partial  current  problem  description  and  implementation  for  both  the  zeroth  and 
first  moment  methods  are  also  described  in  chapter  four.  The  validation  of  the 
code  is  presented  in  chapter  five  as  well  as  testing  designed  to  demonstrate  that 
the  DI  method  performs  at  least  as  well  as  other  methods  for  routine  problems. 

The  DI  method  is  tested  on  a  variety  of  problems  in  the  remaining 
chapters.  The  testing  is  designed  to  illustrate  two  points:  the  DI  method 
performs  at  least  as  well  as  other  methods  and  the  DI  method  works  for  those 
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problems  where  other  methods  fail  or  have  difficulties.  The  testing  is  also 
designed  to  determine  the  limitations  of  the  DI  method.  The  problems  where 
other  methods  have  difficulties  are  presented  in  chapter  six  and  problems  that 
stress  the  DI  method  are  presented  in  chapter  seven. 
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II.  Theory 


A.  The  Discrete  Ordinates  System  of  Equations 

The  discrete  ordinates  system  of  equations  is  derived  from  the  linear  BTE 
by  discretizing  in  space  and  angle.  The  system  of  equations  may  be  expressed  (9: 
166)  as: 

tin  ■  V  y/(r,Q.n )  +  <7t (r)y/(r, Q„)  =  S(r, Q „)  +  E(r, Q„)  ,  (2.1) 

with  an  appropriate  spatial  discretization.  This  results  in  a  system  of 
simultaneous  equations  that  is  too  large  to  solve  directly.  Therefore,  this  system 
of  equations  is  solved  by  source  iteration.  Advances  in  computing  speed  and 
available  memory  suggest  another  approach  is  appropriate,  motivating  this 
research. 

Rather  than  try  to  solve  the  large  problem  directly,  the  intent  is  to  break  a 
single  large  problem  into  two  smaller  problems  that  can  be  coupled  together. 

The  two  problems  can  be  described  as  a  local  detailed  balance  problem  in  each 
cell  of  a  spatial  grid  and  a  global  flow  balance  problem.  Both  problems  assemble 
the  discrete  ordinates  system  of  equations  in  a  form  that  gives  the  outgoing 
particle  flow  in  terms  of  the  inward  particle  flow  and  emissions.  By  determining 
the  proper  balance  on  both  scales,  using  the  local  balance  to  improve  the 
coefficients  in  the  global  balance  equations,  and  the  solution  to  the  improved 
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global  balance  to  improve  the  local  balances,  the  problem  can  be  solved 
iteratively. 

In  general,  a  cell  system  of  equations  is  written: 

r=K0,7"  +KosS  +  KoeE,  (2.2) 

^  =  Kr,T+K^S  +  KrEE,  (2.3) 

and 

S  =  (2.4) 

-^out 

In  equations  (2.2)  through  (2.4),  j  is  a  vector  of  coefficients  of  basis 
functions  in  an  approximation  to  the  distribution  of  current  on  the  faces  of  the 

— in 

cell  for  the  outward  directions,  j  is  a  vector  of  coefficients  of  basis  functions  in 
an  approximation  to  the  distribution  of  current  on  the  faces  of  the  cell  for  the 
inward  directions,  y/  is  a  vector  of  coefficients  of  basis  functions  in  an 
approximation  to  the  distribution  of  the  angular  flux  within  the  cell,  S  is  a 
vector  of  coefficients  of  basis  functions  in  an  approximation  to  the  distribution  of 
the  scattering  source  within  the  cell,  E  is  the  vector  of  coefficients  of  basis 
functions  in  an  approximation  to  the  distribution  of  emissions  within  the  cell. 

The  matrices,  K OI ,  K os  ,  K ()E  ,  K(///  ,  K l//S  ,  and  K tj/E  are  the  relations  between 
the  vectors  for  the  spatial  quadrature.  The  matrix  X  contains  the  scattering 
contribution  and  angular  quadrature  weights  to  relate  the  scattering  source  and 
the  angular  flux.  These  equations  are  developed  further  in  Chapters  three  and 
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four.  The  particle  flow  variable  at  the  cell  faces  in  equations  (2.2)  and  (2.3)  could 
be  expressed  either  as  angular  flux  or  as  angular  current.  The  reason  for  using 
currents  is  presented  in  section  B.  These  equations  are  used  for  the  local  detailed 
balance  problem  within  a  cell. 

The  global  balance  problem  uses  the  flow  of  particles  across  cell  faces  for  all 
the  cells  in  the  mesh,  and  removes  the  angular  dependence  by  integrating  a 
modification  to  equation  (2.2)  over  a  hemisphere  to  determine  the  particle  flow. 

A  discussion  of  both  problems  follows. 

1.  Local  Detailed  Balance 


Local  detailed  balance  is  found  by  eliminating  the  scattering  source  from 
the  system  of  equations  in  a  cell.  This  allows  the  direct  calculation  of  the 
detailed  flow  of  particles  in  a  cell  from  the  flow  from  adjacent  cells  and  emissions 
within  the  cell,  again  accounting  for  all  of  the  scatters  a  particle  can  undergo 
within  the  cell.  The  local  detailed  balance  relation  for  a  cell  is: 


-Out  -In  — 

J  =m  oil  +m  oeE. 


(2.5) 


Again,  j  is  the  current  at  a  cell  edge  for  all  the  ordinates  in  the  angular 


quadrature  set,  E  is  the  emissions  in  the  cell  along  each  ordinate,  and  mOI  and 
m()E  are  matrices  which  give  the  contributions  of  the  inward  particle  flow  and 
emissions  respectively. 


To  convert  equations  (2.2)  through  (2.4)  into  the  form  of  equation  (2.5), 


substitute  equation  (2.4)  into  equation  (2.3)  and  solve  for  the  angular  flux: 
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(2.6) 


—  -*in  — 

^  =  LK  +  LK^F, 

where 

L  =  (I-Kt,sSr1  (2.7) 

This  result  and  equation  (2.4)  are  substituted  into  equation  (2.2)  for  the  current: 

T‘  =KoiT  +KOSnLK¥i7"  +LKveE)  +  KoeE  .  (2.8) 

Collecting  terms  yields: 

-~out  —  in  — 

j  =  (K0£+K0S£LK,,;X/  +(K0£+KosSLK„£)£.  (2.9) 

The  matrices  in  the  parentheses  are  in  the  form  of  equation  (2.5).  Further,  the 
first  matrix  represents  the  contribution  to  first  flight  of  particles,  while  the 
product  term  represents  the  contribution  from  particles  after  scattering.  The 
convention  of  bold  symbols  represents  matrices,  while  lower  case  m  is  a  reminder 
that  this  is  a  cell  formula.  This  provides  the  needed  formulas  for  the  coefficient 
matrices: 

mo;=K0;+KosZLKr/,  (2.10) 

and 

mOE  =  +  KOS  SLK  rE.  (2.11) 

2.  Global  Flow  Balance 

The  global  flow  problem  solves  directly  across  the  problem  for  the  flow  of 
particles  across  cell  edges  with  no  angular  dependence.  Equation  (2.5)  is 
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integrated  over  the  appropriate  angles  to  determine  the  outward  particle  flow  for 


a  cell  and  this  is  used  to  create  a  system  of  equations  across  the  spatial  mesh: 

A  x  =  b,  (2.12) 

where  the  flow  of  particles  across  cell  edges,  x  is  only  dependent  on  the  forcing 
term,  b  . 

The  global  flow  problem  is  much  smaller  than  the  discrete  ordinates 
system  of  equations.  The  matrix  A  is  so  sparse  that  the  system  can  be  solved 
directly.  The  detailed  angular  information  for  the  flow  of  particles  is  implicit  (in 
the  elements  of  A);  the  partial  currents  of  particles  passing  through  all  cell  edges 
are  the  only  (explicit)  unknowns.  This  uses  the  spatial  quadrature  to  model  the 
contribution  of  particles  entering  the  cell  from  any  direction,  scattering  any 
number  of  times,  and  exiting  the  cell  edge. 

The  cell  flow  of  particles  in  the  local  balance  problem  contains  the  detailed 
angular  information  that  is  implicit  in  the  global  balance  solution.  However,  the 
detailed  cell  flow  does  not  necessarily  include  the  contribution  from  particles  that 
flow  from  nonadjacent  cells  after  any  number  of  scatters.  To  overcome  this 
shortcoming  and  retain  angular  information  for  the  global  flow  problem  the  two 
problems  must  be  linked.  Coupling  the  global  and  detailed  balance  problems 
solves  this. 
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3.  Coupling  the  Local  Balances 


The  cell  local  balances  are  coupled  across  the  spatial  mesh  and  with  the 
global  balance  problem  through  cell  coupling.  The  level  of  particles  scattering 
among  cells  is  contained  in  the  global  balance  solution,  but  the  distributions  of 
particles  in  angle  is  not.  Coupling  the  local  balances  addresses  the  role  of  these 
distributions.  This  allows,  through  an  angular  integration,  the  relative 
importance  of  an  angular  direction  to  be  used  in  the  global  balance  solution  and 
the  distribution  of  the  appropriate  level  for  a  particular  direction  back  to  the 
local  detailed  balance  problem  from  the  global  balance  solution. 

If  the  correct  coupling  were  known,  both  problems  could  be  solved  exactly 
and  the  detailed  and  global  flow  of  particles  could  be  calculated  directly.  As  it  is 
not  known,  an  estimate  is  used  and  iteration  is  used  to  improve  the  estimate  of 
the  coupling.  This  is  not  source  iteration;  instead,  this  iteration  seeks  to  improve 
the  estimate  of  the  coupling  of  the  local  balances  rather  than  improving  the 
scattering  source  estimate.  The  general  method  is  shown  in  figure  2.1. 

The  figure  shows  the  general  phases  of  the  distribution  iteration  method. 

An  assumption  is  made  for  cell  edge  distributions.  This  assumption  is  used  to  set 
up  and  solve  the  global  flow  problem,  taking  into  account  emissions  within  the 
problem  and  boundary  conditions  to  set  approximate  cell  edge  flow  values.  The 
cell  edge  values  are  used  as  inflows  for  neighboring  cells  to  find  cell  outflows  and 
improve  the  estimate  of  the  edge  distributions. 
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Better 

Detailed 

Distributions 


Figure  2.1  Flowchart  describing  the  general  distribution  iteration  method. 


This  inflow  to  outflow  can  be  repeated  until  the  estimate  of  the  edge 
distribution  is  sufficiently  improved.  The  innermost  loop  shows  where  additional 
iteration  is  done  for  nonlinear  spatial  methods,  which  are  not  included  in  this 
research.  Better  inward  and  outward  detailed  particle  flow  solutions,  through 
iteration  if  needed,  provide  better  coupling  with  detailed  distributions.  The 
updated  edge  distribution  is  used  to  set  up  another  global  flow  problem  and 
improve  the  estimate  for  the  global  flow  problem  solution,  and  the  process 
repeats  until  a  convergence  criterion  is  met. 
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B.  Choices  that  Define  a  Method 


The  discrete  ordinates  system  of  equations  is  defined  by  the  angular  and 
spatial  discretizations  that  are  used.  These  choices,  along  with  the  possible 
implementation  choices  shown  in  figure  2.1  define  a  distribution  iteration 
method.  A  review  of  the  choices  and  considerations  for  implementation  follow. 

1.  Angular  Quadrature  Sets 

In  general,  the  angular  quadrature  sets  are  used  to  evaluate  the  angular 
integrals  needed  for  solutions  to  the  discrete  ordinate  equations, 
a.  Slab  Geometry  Angular  Quadrature  Sets 

In  slab  geometry,  two  different  quadrature  sets  were  considered.  A  brief 
description  of  these  quadratures  follows. 

Discrete  Elements  Quadrature  Set 

Discrete  elements  (DE)  quadratures  were  used  by  Wager  to  demonstrate 
the  Angle  Iteration  method  in  one  dimension  (16:  2-3-5).  The  discrete  elements 
quadratures  do  not  exactly  compute  the  factor  of  1/3  in  the  diffusion  coefficient. 
The  angular  quadrature  should  exactly  integrate  the  following  integral: 

j^V=b  (2.13) 

-1  1  3 

where  LI  represents  the  direction  cosine.  In  slab  geometry  for  a  discrete  elements 
quadrature  set  with  an  even  number,  N,  of  equal  weight  elements,  the  general 
expression  for  the  mean  jUn  in  an  element  is: 
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(2.14) 


Mn=^~ 


(2/7-1) 

N 


where  the  element  size  is: 


AjU  =  — . 

1  N 


(2.15) 


This  particular  case  for  discrete  elements  is  akin  to  a  composite  mid-point 
method.  Using  the  DE  angular  quadrature,  the  integral  in  equation  (2.13)  is 


approximated  as: 


f  d/U  1  |  r)  1.  1 

“2  HVnAM„=- 

n= 1 


-1 


1 - - 

V  N*J 


3 


(2.16) 


As  shown  in  equation  (2.16),  higher  resolution  DE  quadratures  (larger  N)  are 
closer  to  meeting  the  diffusion  limit,  but  are  still  not  exact. 

While  this  quadrature  set  made  the  visualization  of  collapsing  and 
allocating  the  angular  flux  easier,  the  discrete  elements  quadratures  do  not  meet 
the  diffusion  limit.  Problems  which  are  highly  scattering,  which  are  the  type  of 
problems  where  synthetic  acceleration  and  source  iteration  have  difficulty,  and 
that  this  research  will  examine,  need  an  angular  quadrature  that  meets  the 
diffusion  limit. 

Discrete  Ordinate  Quadrature  Sets 

In  the  case  of  discrete  ordinates  angular  quadratures,  an  exact  relation  for 
the  integral  in  equation  (2.13)  is  often  considered  to  be  a  requirement  for  useful 


quadratures.  Lewis  and  Miller  provide  a  description  of  common  discrete 
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ordinates  angular  quadrature  sets  for  slab  geometry  (9:  119-126)  and  XY- 
geometry.  In  general,  these  are  even  quadratures  which  are  symmetric  about 
fj.  -  0  since  positive  and  negative  particle  flows  are  generally  of  equal  importance. 
A  discussion  of  two  popular  quadrature  sets  for  slab  geometry  follows. 

Single  Range  Gauss-Legendre 

This  quadrature  set  is  also  known  as  PN  quadratures.  The  ordinates, [ln  , 
are  the  N  roots  of  the  Legendre  polynomial: 


^iv(/L)  =  0’  n  =  l,2,...,N. 


(2.17) 


The  weights  are  found  such  that  the  quadrature  set  correctly  integrates  all 
polynomials  through  order  2N  —  1.  The  symmetry  of  the  PN  quadrature  set  and 
the  properties  of  the  Legendre  polynomials  make  this  quadrature  set  popular  for 
certain  problems  (9:  119-121). 

Double  Range  Gauss-Legendre 

Double  range  or  DPN  quadrature  sets  (9:  121-126)  are  similar  to  PN 
except  that  quadratures  are  developed  for  the  integrals  over  - 1  <  //  <  0  and 
0  <  )l  <  1 .  These  quadrature  sets  are  used  for  their  improved  treatment  of 
vacuum  boundaries.  In  our  case,  these  are  desirable  because  partial  currents  are 
defined  as  integrals  over  these  two  domains.  For  example,  consider  the  function 
defined  as: 


f(M)  =  0  -l<ju<0, 

/(//)  =  //  0  <  //  <  1. 


(2.18) 
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The  integral  over  the  whole  domain  is 


\f(n)dn  =  \.  (2.19) 

-l  2 

Double  range  quadrature  sets  can  integrate  this  function  exactly  (even  for  the 
lowest  order),  while  single  range  quadratures  do  not.  Table  2.1  shows  the 
integration  results  for  several  single  range  quadrature  sets  for  equation  (2.19)  and 
demonstrates  this.  The  single  range  quadratures  of  order  2  -12  and  the 
integration  results  of  equation  (2.19)  are  shown.  The  error  listed  in  the  last- 
column  is  the  absolute  difference  from  the  exact  solution. 


Table  2.1  Single  range  Gauss  quadrature  results  for  equation  (2.19) 


Pn 

Integration 

Results 

Error 

(Difference) 

Pi 

0.57735 

0.07735 

Pa 

0.52126 

0.02126 

Pe 

0.50994 

0.00994 

P& 

0.50576 

0.00576 

Pio 

0.50376 

0.00376 

P\2 

0.50264 

0.00264 

b.  XY  -  Geometry  Angular  Quadrature  Sets 

The  discrete  ordinates  quadratures  were  implemented  to  evaluate  the 
angular  integrations  needed  for  solutions  to  the  discrete  ordinates  equations  in 
XY  -  geometry  because  they  meet  the  diffusion  limit.  Two  different  quadratures 
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were  implemented  in  order  to  compare  with  previously  published  results  and  to 


show  the  facility  of  the  method  in  implementing  different  angular  quadratures. 

Level  Symmetric 

Level  symmetric  quadratures  are  widely  used  (9:  158-162),  and  this 
quadrature  set  was  used  to  compare  with  published  results.  These  quadrature 
sets  are  referred  to  as  SN  quadratures  and  contain  the  same  set  of  direction 


cosines  with  respect  to  each  axis. 


There  are 


N(N  +  2) 
8 


ordinates  per  octant. 


The 


quadrature  weights  meet  the  condition  that  all  weights  must  be  equal  for  points 
obtained  by  permuting  the  direction  cosines.  A  useful  property  of  the  SN 
quadratures  is  that  the  ordinate  directions  are  invariant  to  90°  rotations  about 
any  axis.  The  quadratures  sets  S4,  S6,  and  Ss  were  implemented  for  XY  - 
geometry. 


Product  Quadratures 

A  product  quadrature  was  also  implemented  to  show  the  facility  with 
which  different  angular  quadratures  could  be  used.  Abu-Shumays  (1:  299-301) 
showed  that  a  quadruple  range  quadrature  set  was  competitive  for  improving 
accuracy.  In  this  quadrature  method,  the  polar  angle,  </> ,  is  integrated  using  a 
Gauss- Cristoffel  quadrature  and  the  azimuthal  angle,  (0 ,  is  integrated  using  a 
Gauss-Chebychev  quadrature.  The  direction  cosines  are  calculated  using  these 
two  quadratures: 
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nn  =cos(«m(lj))sin(^(>j)), 

%=sin(«m(„))sin(^(„)), 

and  the  final  quadrature  weights  are  the  product  of  the  quadrature  weights, 

wn  =  wm(n)wf(n)  > 


(2.20) 


(2.21) 


for  both  the  Chebychev  and  Cristoffel  quadratures.  Cristoffel  quadratures  with 
1-3  levels  per  octant  and  Chebychev  quadratures  with  2-5  levels  per  octant  were 
implemented  for  XY  -  geometry. 

2.  Spatial  Quadratures 


Spatial  quadratures  methods  can  be  characterized  by  the  highest-order 
spatial  moment  balance  that  is  satisfied  exactly.  Zeroth-moment  and  first- 
moment  methods  are  used  here.  An  advantage  of  linear  methods  is  that  for  the 
distribution  iteration  methods,  the  matrix  relationships  providing  the  flow  of 
particles  and  defined  by  the  angular  and  spatial  quadratures  are  fixed  and  do  not 
need  to  be  calculated  for  each  iteration.  For  this  reason,  only  linear  methods  are 
used  here. 

Several  attributes  of  spatial  quadrature  methods  are  of  especial  interest: 
positivity,  linearity,  and  (2nd  order  or  better)  accuracy.  Positivity  means  that  the 
outgoing  face  flow  value  (and  the  flux  within  the  cell)  returned  by  the  spatial 
method  is  nonnegat-ive,  given  nonnegative  inflow  flux  and  source.  Negative  flow 
values  are  non-physical  and  are  strictly  an  artifact  of  the  spatial  method. 
Linearity  refers  to  the  superposition  of  solutions,  a  solution  for  a  source  that  is 
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the  sum  of  other  sources  is  also  the  sum  of  the  solutions  for  the  other  sources  (11: 
33).  Accuracy  refers  to  the  truncation  error  on  fine  meshes  (9:  371).  A  spatial 
method  has  at  most  two  of  the  three  attributes  (9:  135).  The  description  of  how 
these  attributes  align  with  the  choice  of  spatial  quadrature  follows. 

Step  characteristic  (SC)  is  a  zeroth  spatial  moment  method.  It  is  a  linear 
and  positive  method,  but  it  is  1st  order  accurate.  The  SC  method  was 
implemented  for  slab  and  XY-geometry  to  demonstrate  the  method  and  to 
compare  with  Wager’s  results. 

Weighted  diamond  difference  (WDD)  is  a  zeroth  spatial  moment  method 
that  is  also  a  linear  and  positive  method,  but  has  less  than  2nd  order  accuracy. 

The  method  is  used  in  production  codes  and  was  used  by  Azmy  to  demonstrate 
the  loss  of  effectiveness  for  DSA.  The  WDD  method  was  implemented  in  XY- 
geometry  to  compare  with  published  results. 

Linear  discontinuous  (LD)  is  a  first  spatial  moment  method  that  is  linear 
and  3rd  order  accurate,  but  is  not  a  positive  method.  The  LD  method  is  also  used 
in  production  codes  and  is  one  of  the  spatial  quadratures  that  meets  the  diffusion 
limit  on  thick  cells.  The  LD  method  was  implemented  for  slab  and  XY-geometry 
to  demonstrate  the  DI  method  for  first  spatial  moment  methods. 

Linear  characteristic  (LC)  is  another  non-positive  first  spatial  moment 
method  that  is  linear  and  4th  order  accurate.  It  is  used  for  better  accuracy,  but 
does  not  (like  all  characteristic  methods)  meet  the  thick-cell  diffusion  limit.  The 
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LC  method  was  implemented  for  XY-geometry  to  demonstrate  the  flexibility  of 


the  DI  method  for  first  spatial  moment  methods. 

b.  Nonlinear  Methods 

Wager  demonstrated  the  feasibility  of  using  the  exponential  characteristic 
(EC)  method  for  his  work  in  slab  geometry,  a  nonlinear  method  (16:  6-1). 
Nonlinear  methods  add  additional  complexity,  but  have  the  attributes  of 
accuracy  (EC  has  4th  order  accuracy),  ability  to  use  a  coarser  spatial  mesh  for 
accurate  solutions,  and  positivity.  However,  additional  calculations  are  needed 
for  the  innermost  loop,  as  noted  in  Figure  2.1.  Due  to  the  additional  complexity 
required  for  nonlinear  methods,  implementation  of  DI  with  EC  is  left  for  future 
efforts. 

3.  Cell  Face  Flow  Variables 

Two  choices  for  cell  face  flow  variable  are  readily  apparent:  angular  flux, 
!// .  and  angular  current,  j  . 

a.  Angular  Flux  as  Cell  Face  Variable 

Angular  flux  is  commonly  used  for  the  cell  face  flow  variable.  Spatial 
quadratures  are  presented  in  the  literature  in  terms  of  angular  fluxes.  However, 
with  this  choice,  the  global  flow  balance  variable  lacks  physical  meaning;  it  is  the 
angular  integral  of  the  angular  flux  over  a  hemisphere,  which  is  neither  a  partial 
current  nor  a  scalar  flux. 
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b.  Angular  Current  as  Cell  Face  Variable 


Angular  current  has  not  been  used  for  the  cell  face  flow  variable  (to  my 
knowledge),  but  can  be  determined  easily  from  the  angular  flux  and  the  direction 
cosines  from  the  angular  quadrature.  The  angular  current  is 

in /out  ^  *  in  I  out  ^ 

j face  ~  ^ '  n  face  ^(Q)  ,  where  in  or  out  is  chosen  for  a  given  £1  such  that  the  dot 

product  is  positive.  Also,  the  global  flow  variable,  ,  is  now  a  physical  quantity: 
the  partial  current  through  the  cell  face.  This  changes  the  global  flow  problem  to 
a  partial  current  problem  that  is  an  explicit  statement  of  conservation  of  particles 
within  each  cell.  This  motivated  my  choice  of  face  flow  variable  for  the 
distribution  iteration  method.  The  difference  is  more  than  one  of  bookkeeping; 
distribution  iteration  using  the  partial  current  problem  converges  in  fewer 
iterations,  as  demonstrated  by  the  testing  presented  in  chapter  3. 

4.  Coupling  the  Local  Balances 

Coupling  the  local  balances  requires  carrying  information  about  particle 
flow  from  each  cell  in  the  spatial  mesh  to  other  cells,  including  cells  that  are  not 
adjacent  and  may  be  distant.  In  order  to  obtain  a  rapidly  converging  method,  an 
efficient  coupling  method  is  needed.  Three  different  options  are  presented:  local 
balance  sweeping;  red/black;  and  discrete  ordinates  sweeping, 
a.  Local  Balance  Sweeping 

In  the  local  balance  sweeping  method,  the  current  cell  uses  the  outflows 
from  the  adjacent  cell  as  inflows.  This  method  is  easy  to  implement,  but  requires 
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multiple  sweeps  to  communicate  between  nonadjacent  cells.  Additionally,  the 


outflow  in  the  direction  of  the  sweep,  which  is  the  inflow  for  the  next  cell  has 
more  improvement  than  the  inflow  at  far  side  of  the  next  cell.  There  is  a 
possibility  that  the  inflow  estimate  on  the  cell  edges  may  not  have  equal 
improvement  across  the  cell,  which  may  introduce  a  bias  across  the  cell.  This  is 
the  method  Wager  used  in  his  efforts  (16:  2-68).  The  sweep  is  sequential  in  Id, 
and  has  some  parallelism  in  higher  dimensions, 
b.  Red/Black 

The  red/black  method  divides  the  spatial  mesh  into  alternating  cells  and 
assigns  a  color,  similar  to  a  checkerboard.  All  the  red  cells  can  be  done  in 
parallel.  The  red  cell  outflows  are  the  black  cell  inflows  so  that  the  black  cells  can 
then  be  done  in  parallel.  Each  cell  communicates  only  locally  -  to  its  immediate 
neighbors  in  the  spatial  grid.  This  is  the  ideal  situation  for  fully  parallel 
computations  with  efficient  scaling  to  many-processor  systems. 

However,  the  region  influenced  by  a  localized  source  in  a  problem  with 
little  scattering  is  extended  by  only  one  cell  (in  all  directions)  for  each  red  or 
black  calculation,  hence  two  cells  per  red/black  iteration.  Thus  convergence  may 
be  slow  for  such  cases.  These  are  the  conditions  in  which  SI  works  best,  because 
the  sweeps  along  the  ordinate  carry  the  first-flight  influence  of  a  localized  source 
throughout  the  problem  in  one  iteration.  This  motivated  the  next  approach. 
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c.  Discrete  Ordinate  Sweeping 


Rather  than  calculate  the  outflow  values  directly,  the  inflow  values  and 
emissions  in  a  cell  can  be  used  to  improve  the  current  estimate  of  the  cell 
scattering  source  including  all  numbers  of  scatters  within  the  cell.  This  way  of 
calculating  the  scattering  source  overcomes  the  difficulties  of  traditional  discrete 
ordinates  methods  of  estimating  the  scattering  source,  particularly  for  high 
scattering  ratios.  The  cell  scattering  source  calculations  can  be  done  in  any 
order,  and  are  also  parallelizable.  The  cell  scattering  sources  are  then  used  for  a 
single  discrete  ordinates  sweep  for  each  ordinate  to  determine  the  cell  outflow 
values.  For  code  implementation,  two  different  discrete  ordinates  sweeping 
methods  were  used.  The  first  was  a  single  source  calculation  followed  by  one 
discrete  ordinates  sweep.  This  proved  sufficient  for  most  problems.  The  other 
method  was  an  adaptive  technique  which  varied  between  one  and  ten  sweeps 
depending  on  the  properties  of  the  problem.  For  each  sweep,  the  scattering 
source  was  updated  using  the  cell  edge  values  and  the  scattering  source  was  used 
to  calculate  new  cell  edge  values.  This  was  used  for  slab  geometry  problems  and 
some  of  the  XY-geomet-ry  problems.  A  further  analysis  and  description  is 
presented  in  chapter  seven. 

A  strength  of  the  discrete  ordinates  sweep  is  that  it  rapidly  communicates 
cell  information  across  the  spatial  mesh  as  the  angular  flux  calculations  in  an 
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ordinate  are  done  over  the  spatial  mesh.  The  sweep  is  limited  in  parallelism,  but 


the  method  may  have  merit  for  serial  machine  implementation. 

d.  Parallel  Efficiency  of  Red/Black  vs.  Sweeping 

Any  algorithm  that  sweeps  through  a  regular,  orthogonal  grid,  such  as  a 
checkerboard  or  its  extension  to  3d,  in  a  compound  direction,  for  example 
upward  and  to  the  right,  is  constrained  by  data  dependencies.  For  example,  after 
the  bottom-left  cell  is  done,  the  data  for  both  the  cell  to  its  right  and  the  cell 
above  it  are  available.  These  two  cells  can  be  done  in  parallel,  after  which  the 
three  cells  above  and/or  right  of  them  can  be  done  in  parallel  and  so  on.  Thus, 
one  sweeps  a  diagonal  line  of  cells  (crosswise  to  the  flow)  in  XY-geomet.ry,  or  a 
diagonal  plane  of  cells  in  XYZ-geometry.  This  is  partially  parallel,  but  much  less 
efficient  than  red/black  (per  iteration).  Let  d  be  the  number  of  spatial 
dimensions  and  n  be  the  size  of  the  mesh  (in  each  dimension). 


Table  2.2.  Parallel  Efficiency  considerations. 


d 

Sweeps 

Stages  per 

Sweep 

Stages 

Asymptotic  Stage  Ratio, 

Sweeps  :  Red/Black 

1 

2 

n 

2/7 

77 

2 

4 

2/t  —  l 

8/7-4 

4/7 

3 

6 

3/7-2 

18/7-12 

9/7 

For  large  n  ,  the  asymptotic  ratio  (large  n  )  of  the  number  of  parallel 
stages  per  iteration  for  sweeping  to  the  number  of  parallel  stages  for  red/black 

2  i 

(two  stages)  per  iteration  is  d  n .  This  analysis  applies  to  both  discrete  ordinates 
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sweeping  and  to  local  balance  sweeping  because  all  the  ordinates  that-  have  the 
same  data  dependency  (such  as  upward  and  to  the  right)  can  be  done  in  parallel 
in  the  diagonal  sweep. 

Consequently,  it  is  reasonable  to  expect  that  red/black  will  be  more 
efficient  for  large  problems  on  MMP  systems  because  the  number  of  red/black 
iterations  should  be  much  less  than  the  grid  size  n  . 
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III.  Slab  Geometry  Implementation  and  Testing 


A.  Local  Detailed  Balance  Problem 

1.  Zeroth  Spatial  Moment  Methods 

Wager  presented  the  foundation  for  the  zeroth  spatial  moment  methods 
(16:  2-22-27)  using  the  angular  flux  formulation.  The  analogous  angular  current 
formulation  is  presented  here.  For  the  zeroth  spatial  moment  methods,  the 
system  of  equations  for  a  cell  are: 


-"■out  — in  a  —  a 

J  =ko/7  +  ^OEAE  5 

(3.1) 

¥A  =  Ku/  j  +  kasa$A  +  kaeaEA  > 

(3.2) 

and 

SA=HSVA- 

(3.3) 

In  these  equations,  K0/  ,  K OSA  ,  K OEA  >  K A[  ,  K ASA  and  K  AEA  ,  are  diagonal 
matrices  of  transport  coefficients  that  define  the  relations  of  the  inputs  of  a  cell 
to  the  calculated  quantity.  For  example,  K0/  represents  the  contribution  to  the 
outgoing  flux  from  the  incoming  flux  and  K ASA  represents  the  contribution  to 
the  average  flux  from  the  average  scatter.  The  values  of  the  transport 
coefficients  are  determined  from  the  spatial  quadrature  used.  Letting  D(x)  be 
the  diagonalization  operator  that  creates  a  diagonal  matrix  from  vector  x ,  the 
general  matrices  become: 

K  =  D(ifc)  .  (3.4) 
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-*out  — in  — -A 

The  quantities  j  ,  j  and  !//  are  the  variables  to  represent  the  cell  outward 
angular  current,  inward  angular  current  and  average  angular  flux  respectively. 
The  zeroth  spatial  moment  over  a  cell  is  normalized  to  be  the  average  value. 
Thus,  the  quantities  SA  and  EA  are  the  variables  to  represent  the  average 
scatter  and  average  emissions  in  a  cell.  The  vector  notation  represents  an  array 
for  the  variable  with  all  the  ordinates  in  the  angular  quadrature.  Also,  for 
isotropic  scatter: 

Z5=<ts1D(w),  (3.5) 

where  D(w)  is  the  diagonal  operator  on  quadrature  weight  vector,  <JS  is  the 
isotropic  scattering  source  and  1  is  a  matrix  with  one  for  every  element.  In 
general,  this  matrix,  Z5  contains  the  scattering  cross  sections  with  the 
appropriate  quadrature  weights  to  calculate  the  scattering  source  from  the 
average  flux. 

Equations  (3.3)  can  be  substituted  into  equation  (3.2)  to  solve  for  the 
average  flux  in  a  cell  in  terms  of  the  incoming  cell  angular  current  and  the 
emission  in  a  cell: 

¥A  =(I-K^Zsr1K^r+(I-KJMZs)"1KJ£J£".  (3.6) 

Equations  (3.6)  and  (3.3)  can  be  substituted  into  equation  (3.1)  to  solve  for  the 
outgoing  angular  current  in  a  cell  in  terms  of  the  incoming  cell  angular  current 
and  the  emission  in  a  cell: 
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(3.7) 


j  ~  (Ko/  +K0&4  Zs(I ~KAsa  1^Ai)j  + 

(K0£4  +K0&4  £s(I-KA£/4  ^5)  1k^£^)^j4- 

Again,  the  factor  (I-K^^Xs)-1  can  be  thought  of  as  modeling  infinite  within 
cell  scattering  (16:  2-55).  The  terms  in  the  outer  parentheses  can  be  expressed  as 
a  single  matrix: 

-out  -in  —  a  .  . 

j  =  moiJ  +  ™oee  ■  (3-8) 

An  exactly  analogous  derivation  is  done  for  the  angular  flux  formulation.  For  the 
zeroth  spatial  moment  methods,  these  matrices  need  only  be  calculated  one  time 
for  each  material  (with  a  uniform  spatial  mesh).  Equation  (3.8)  can  be  used  to 
solve  the  cell  detailed  balance  problem  for  both  the  local  balance  sweeping  and 
the  red/black  methods.  Equations  (3.6)  and  (3.3)  are  used  to  determine  the  cell 
scattering  sources  for  the  discrete  ordinates  sweep  method. 

Zeroth  Moment  Transport  Coefficients 

For  the  testing  in  slab  geometry,  the  transport  coefficients  for  the  step 
characteristic  will  be  discussed  for  zeroth  spatial  moment  methods.  The 
transport  coefficients  are  used  to  build  the  diagonal  matrices  used  in  the  local 
detailed  balance  problem.  Other  zeroth  spatial  moment  methods  would  be 
implemented  using  the  same  procedure. 
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1.  Angular  Flux  Formulation 


Lathrop  introduced  the  step  characteristic  method  in  1969  (11:  24)  and 
the  quadrature  equations  can  also  be  found  in  the  literature  (10:  v-8).  Wager 
presented  SC  relations  (16:  2-18,2-22)  in  a  more  compact  notation  using  the 
exponential  moment  functions  of  order  m  developed  by  Mathews  et  al.  (11:  27) 
where: 

l 

Mm(x)  =  jdt(l-t)me~xt.  (3.9) 

o 


The  cell  optical  thickness  measured  along  ordinate  n  is  used  in  these  relations 
and  is  defined  as: 


oAx 


(3.10) 


where  <7  is  the  total  cross  section,  Ax  is  the  cell  width  and  jUn  is  the  direction 
cosine  from  the  angular  quadrature.  As  an  example,  consider  a  quadrature  set 
with  four  ordinates,  (1  and  2  to  the  right,  3  and  4  to  the  left).  The  cell  SC 
equations  for  the  outgoing  angular  fluxes  are: 

Vi=e~£lVi 

l/f?  =  0~£2  1//T  - 


(3.11) 

W  N 

(3.12) 

A-M0U,)S,a  +—M0(e,)E? , 

\Mi\  Mi, 

(3.13) 
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and 


¥a  =  e~£^l  +^-,M0{£A)St  +  —M0{e,)Et , 
N  Ea 


(3.14) 


where  R  and  L  designate  the  right  and  left  faces  of  the  cell.  Equations  (3.11) 
through  (3.14)  take  the  form  of  equation  (3.1),  with  K0/  =D (koi) , 


K 


OSA  =D (koSA)  and  K Qea  =D (/coea)  ,  where  the  k  vectors  are: 


and 


(^0/)«  e  "  1 

(3.15) 

Ax 

(W,=nMo(f»)> 

(3.16) 

\e„\ 

Ay 

(kOEA\=—M0(£n). 

(3.17) 

En 

For  the  same  quadrature  set,  the  SC  equations  for  the  average  angular  flux  are: 


Wi  =  mq(£\)¥\L  +—Ml(£l)SlA+—Ml(£l)ElA  ,  (3.18) 

Ei  Ei 


¥i  =  M0(£2)¥2  +—M](£2)Sa  +—M1(£2)Ea  ,  (3.19) 

Ei  Ei 


y/3  =  M0(£3)y/3  +—Ml(£3)SA+—Ml(£3)E^  ,  (3.20) 

Ei  Ei 


and  yi  =  M0(£4)itf  +-^Ml(£4)St  +—M1(£4)e£  .  (3.21) 

\Ea\  Ea 


Equations  (3.18)  through  (3.21)  take  the  form  of  equation  (3.2),  with 

Kai  -  D(k/i;  ) ,  K ASA  =D (kAsA)  and  KAEA  =  D(k  aea  ) ,  where  the  k  vectors  are: 

Vai\=M: 0(e„),  (3.22) 
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Ax 

(k ASA )n  -  i  I  M\{en)  , 

(3.23) 

\Mn\ 

and 

Ax 

(k aea)h - 

(3.24) 

Mn 

2.  Angular  Current  Formulation 

The  x  component  of  current  along  an  ordinate  in  slab  geometry  is  jn  —  \m,\  ¥n  ■ 

The  general  form  for  the  angular  flux  equations  is  exactly  analogous  to  the 

angular  current  equations  (3.1)  through  (3.3),  but  the  formulas  for  some  of  the  k 

vectors  are  different.  To  change  the  transport  coefficients  to  the  current 

representation  requires  multiplying  or  dividing  the  transport  coefficient  by  |//„| 

where  appropriate.  The  system  of  equations  for  the  outgoing  flow  and  cell 

average  flux  shown  in  equations  (3.11)  through  (3.14)  and  equations  (3.18) 

through  (3.21)  must  also  be  changed.  Multiply  equation  (3.11)  by  |//,  to  get  the 

corresponding  angular  current  formulation  equation: 

jf  =  e~£'AL  +  A xMQ(ex)Sf  +  AxMoOq^f4 .  (3.25) 

Thus  the  k  vectors  for  the  outward  currents  are: 

(k0l)n=e-£»,  (3.26) 
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(k(JSA  )n  -  AxM0  (£•„)  , 


(3.27) 


and 


(^OEa)h  -  AxM0(en)  . 


(3.28) 


■L  / 

Replacing  y/\  with  the  equivalent  ^  /  i  in  equation  (3.18)  yields  the 

/  I A  I 

corresponding  angular  current  formulation  equation: 


A 


A 


A 


(3.29) 


hence,  the  k  vectors  for  the  average  flux  are: 


and 


(^ AlXi  |  | 

|a«| 

(3.30) 

Av 

(k ASA )n  ~  1  1  M\(£n)  , 

(3.31) 

|Ai| 

At 

(kAEA)n=—Mx(en). 

(3.32) 

Mn 

2.  First  Spatial  Moment  Methods 


The  slab  geometry  zeroth-moment  methods  are  in  the  general  form  of 
equations  (2.3)  through  (2.5);  the  only  difference  is  the  addition  of  “A”  to  some 
of  the  subscripts  and  superscripts.  The  first-moment  methods  are  also  of  that 
form,  as  is  shown  in  this  section.  First,  however,  the  first-moment  methods  are 
presented  in  a  form  that  follows  naturally  from  the  way  such  spatial  quadratures 
are  normally  written.  This  requires  several  additional  terms  and  equations  to 
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account  for  the  contribution  to  the  flux  from  the  first  moment  of  the  scattering 


source: 


j  ~  Ko/  J  +  K<jsaSA  +  K oeaEA  +  KQexEX  > 

(3.33) 

VA  =K4/7"  +KasaSA  +KasxSX  +K aeaEa  +K aexEx  , 

(3.34) 

YX  =  Kx/  j  +  KxsaSa  +  K xsxSx  +  +  K xeXEx  , 

(3.35) 

SA=ZSVA, 

(3.36) 

and 

SX  =  Ys¥X- 

(3.37) 

The  new  matrices,  K OEX  ,  K AEX  •  K Xi  •  K XS4  ,  K XEA  ,K Xsx  and  ^xex  ■ 

are  also 

diagonal  matrices  of  transport  coefficients  that  define  the  relations  of  the  inputs 
of  particles  to  the  calculated  quantity.  For  example  KX!  represents  the 
contribution  to  the  first  spatial  moment  of  the  flux  in  the  cell  from  the  incoming 
flux.  The  diagonal  values,  or  transport  coefficients,  are  also  determined  from  the 

spatial  quadrature  used.  The  new  quantities  !//  ,  SA  and  E  are  the  variables 
to  represent  the  x-moments  of  the  angular  flux,  scatter  and  emissions 
respectively. 

This  system  takes  the  general  form  by  collecting  the  zeroth  and  first 
moments  together  into  single,  larger  vectors.  The  easiest  way  to  do  this  is 
blockwise: 
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¥  = 


¥ 


¥ 


-x 


S  = 


(3.38) 


(3.39) 


and 


(3.40) 


Equation  (3.33)  can  be  written  in  the  general  form 


— out  „  —in  __  —  - 

j  =  K0/  j  +KosS  +  KoeE  , 


by  joining  K  matrices  blockwise: 

K-os  =  [K0&4  K05Z] , 

and  Koe  =  [K  oea  k oex  ]  • 

Similarly,  equations  (3.34)  and  (3.35)  are 

_  —in  —  - 

¥  =  j  +  +  K ¥EE . 


where 


and 


K, 


K„5  = 


K 


.4/ 


1 

II 

x,\’ 

^  ASA 

K  ASX 

J^XSA 

^XSX 

Kara 

^  AFX 

J^XEA 

^  XEX 

(3.41) 

(3.42) 

(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 
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Finally,  equations  (3.36)  and  (3.37)  can  be  also  represented  as: 


S  =  Zys,  (3.48) 

[Ze  0  “ 

where  X  =  (3.49) 

0  Xs 


With  these  matrices  that  include  first  moments,  the  first-moment  system  of 
equations, 


— out  —in  „  —  __  - 

j  =  K0/  j  +  Kos S  +  K oeE  , 

(3.50) 

_  -*Jn  -*  - 

¥~  ^y/I  j  "*"I i 

(3.51) 

and 

031 

II 

M 

(3.52) 

is  of  the  same  form  as  equations  (2.2)  through  (2.4).  Equation  (3.52)  can  be 
substituted  into  equation  (3.51)  to  eliminate  the  scatter: 

¥ = a -iv-  irW + a  -iv  £r‘KV  •  (3-53) 


This  result  with  equation  (3.52)  can  be  substituted  into  equation  (3.50)  to  again 
eliminate  scatter: 


r=(K0;  +Kos  Z(I-K„S  Ir‘Kr/)7"  + 

(Koe  +  Kos  S(I  -  K„s  E)-1  K„£  )E. 


(3.54) 


The  quantities  in  the  parenthesis  can  be  combined  to  form  a  single  matrix 
yielding  the  general  form, 
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(3.55) 


— out  -~in  — 

]  =mOlJ  +  mOEE, 

of  equation  (2.5).  As  with  the  zeroth  spatial  moment  methods,  these  matrices 
are  only  calculated  once  for  each  material  for  first  spatial  moment  methods. 
Again,  equation  (3.55)  can  be  used  to  solve  the  cell  detailed  balance  problem  for 
both  the  local  balance  sweeping  and  the  Red/Black  methods.  Equation  (3.53)  is 
used  to  determine  the  cell  scattering  sources  for  the  discrete  ordinates  sweep 
method.  The  scope  of  this  research  is  limited  to  single  energy  group  problems, 
therefore  the  x-moments  of  emissions  are  zero  and  the  related  transport 
coefficients,  ,K AEX  and  KXEX  ,  were  not  used.  The  code  implementation 

used  an  equivalent  elimination  of  the  scattering  sources,  in  which  a  block  forward 
elimination  and  back  substitution  produces  the  inverse  matrix  in  equation  (3.53). 
Details  are  presented  in  Appendix  E. 

First  Moment  Transport  Coefficients 

For  the  testing  in  slab  geometry,  the  transport  coefficients  for  linear 
discontinuous  methods  are  used  as  an  example  of  a  first-moment,  linear  spatial 
quadrature.  The  transport  coefficients  are  used  to  build  the  diagonal  matrices 
used  in  the  local  detailed  balance  problem.  Other  first  spatial  moment  methods 
would  be  implemented  using  the  same  procedure. 
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Angular  Flux  Formulation 


The  LD  method  has  also  been  used  for  many  years,  and  the  equations  are 
presented  by  Larson  (8:  222)  and  also  by  Lewis  and  Miller  (9:  134).  Similarly  to 
the  zeroth  moment  methods,  a  system  of  equations  for  a  cell  is  set  up  for  all  the 
ordinates  in  the  angular  quadrature  set  for  the  outgoing  flux  in  a  cell,  the  average 
angular  flux  in  a  cell,  and  the  x-moment  of  angular  flux  in  a  cell.  As  was  done 
for  SC,  the  LD  angular  flux  transport  coefficients  are  found  from  the  angular  flux 
equations  in  a  cell,  which  are  shown  in  appendix  A.  The  current  equations  are 
also  derived  in  appendix  A.  The  LD  k  vectors  for  the  outgoing  fluxes  are: 


(koi)n 


6  -  2e„ 


6  +  4  en  +  en 


(3.56) 


(k()SA  h. 


Ax(6  +  en) 


(6  +  4 £„  +£n)  fJ.n 


(3.57) 


C k()Sx)n  ~  77  .  ”  ? 

(6  +  4£„  +  £n  )  //„ 


(3.58) 


and 


(^oea)i  i 


Ax(6  +  £n ) 


(6  +  4£-;i  +£n  )\Mn 


(3.59) 


The  LD  k  vectors  for  the  average  flux  are: 


(k  1  -  6  +  £ n 

\KAl)n  ~ 


6  +  4£n+£n- 


(3.60) 


ASaXi  ~ 


Ax(3  +  £n) 


(6  +  4sn  +£n2)\jun\ 


(3.61) 
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and 


(^ASx)n 


-Ax 

(6  +  4£„  +£n2)\/Jn 


(3.62) 


(^AEa)h  - 


Ax(3  +  £n  ) 

(6  +  4£-„  +£n2)\  fin 


(3.63) 


The  LD  x-moment  angular  flux  k  vectors  are: 


and 


( kxi)n 


(k xsaX ; 


(kxsx)n  ~ 


( kxEA  )r 


_  "3  £„ 

6  +  4en  +  £,~ 

(3.64) 

3Ax 

(3.65) 

(6  +  4en  +£n2)\  fi„\ 

Ax(£n  + 1) 

r,  (3.66) 

(6  +  4  £n+£n2)\fln 

3Ax 

-.  (3.67) 

(6  +  4£n  +£2)\/Jn 

Angular  Current  Formulation 

The  translation  to  the  current  representation  is  done  in  the  same  way  as 
for  the  zeroth  spatial  moment  method  shown  in  equation  (3.25)  using  the  relation 
jn  =  \jUn\l/fn  .  The  LD  k  vectors  for  the  outgoing  currents  are: 


(koi ) 


6-2e„ 


OUn  ,  .  9  5 

6  +  4£-„  +£n 


(3.68) 


(^OSa  )ii  ~ 


Ax(6  +  £n) 
6  +  4£n  +  £n~ 


(3.69) 
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( kosx ) 


Axe„ 


and 


n  s'  a  2 

6  +  4en  +  en 


tu  v  _  Ax(6  +  £„) 
bj  -  7  -  2 

6  +  4f;i  +  £•„ 


The  LD  k  vectors  for  the  average  flux  are: 


(kAi)> 


(k  ASA  )n  ~ 


6 +  £„ 


and 


(k AE/dr 


+ 

VO 

1 

+  £n)\Vn[ 

Ax(3  +  £n ) 

to= 

+ 

VO 

1 

+  £n2)\Mn\ 

- 

-Ax 

to 

+ 

VO 

1 

W)\Mn\ 

Ax(3  +  £n ) 

(6  +  4£n+£n  )\jUn 
The  LD  x-moment-  angular  flux  k  vectors  are: 


(*«)„  =  ic" 


(k xsa)h  ~ 


(kxsx )« 


and 


( kxEA)n  ~ 


(6  +  4£n  +£n  )\/Jn\ 
3Ax 

(6  +  4£n+£n2)\fln\’ 

Ax(£n  + 1) 

(6  +  4£n  +  £n2)\/Jn\ 

3Ax 


(6  +  4£n+£n)\/l„\ 


(3.70) 

(3.71) 

(3.72) 

(3.73) 

(3.74) 

(3.75) 

(3.76) 

(3.77) 

(3.78) 

(3.79) 
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B.  Global  Flow  Balance  Problem 

The  global  flow  balance  problem  determines  the  proper  level  of  flow  values 
across  the  problem.  Setting  up  the  global  flow  balance  problem  with  the  updated 
cell  valued  information  uses  an  array  similar  to  the  flux  weights  used  by  Wager 
(16:  2-73).  However,  a  brief  discussion  of  angular  quadrature  weights  is  needed 
first.  In  the  transport  equations,  the  angular  quadrature  weights  are  used  to 
calculate  the  scalar  flux  and  partial  currents.  In  slab  geometry  the  scalar  flux, 

</>,  is: 

1  1  1 

</>  =  ^\HM)djU~-JZWnV'(lUn)-  (3-80) 

2  _1  2  Vn 

The  partial  currents,  J±  ,  in  slab  geometry  are: 

l 

J+  =  X  WnVnVWi  (3-81) 

0  n3ju„>  0 

0 

and  X  Wn\^’  ’  (3-82) 

-1  n3jU„<0 

where  the  quadrature  weights  are  renormalized  for  each  direction: 

wt  = - ^ - .  (3.83) 

X,  Wn' 
n'3signn’=± 

If  there  are  no  ordinates  where  fln  =  0  ,  (which  is  standard  practice)  and: 

Z  =  1 .  (3-84) 

H3jUn>  0 
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and 


(3.85) 


Z  Wn  =  1  > 

n3jun<  0 

then  =  2,  (3.86) 

Vn 

and  the  quadrature  rules  in  equations  (3.80)  through  (3.82)  are  exact  for 
^(//)  =  constant  (-1  <  //  <  1) .  It  is  sufficient  for  ^  w;;  =  2 ,  no  jtln  =  0  and  a 

\tn 

symmetric  angular  quadrature  set  (which  is  the  case  for  the  angular  quadratures 
tested  in  this  research).  For  these  quadrature  sets,  the  renormalized  weights  in 
equation  (3.83)  are  equal,  wn  =  wn  and  the  wn  notation  is  dropped  for 
convenience. 

The  edge  distribution,  £  ,  is  a  weight  indicating  the  relative  importance  of 
the  current  along  an  ordinate  to  the  partial  current.  The  edge  distribution,  £  ,  is 
defined  for  the  angular  flux  and  current  as  follows: 

C= - - ,  (3.87) 

n  3  Sign  ( //„<  )=Sign  ( ju„  ) 

and  Cl  = - - ,  (3.88) 

Z  Wn'Jn 

n3Sign(fl„-)=Sign(/ln) 

where  wn>  is  the  angular  quadrature  weight.  The  use  of  the  edge  distributions  to 
set  up  the  global  problem  begins  by  referring  to  the  general  form  for  the  cell 
system  of  equations  in  the  current  formulation  in  equation  (3.8): 

-~out  -in  - 

j  =moiJ  +mOEE ■  (3-89) 
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The  denominator  of  equation  (3.88)  is  a  partial  current  as  defined  in  equations 
(3.81)  and  (3.82).  Arranging  the  edge  distributions  and  angular  currents  as 
vectors  for  the  edge  of  the  cell  in  equation  (3.88),  the  inward  angular  current  for 
a  given  direction  is: 

-*in  ~zin  .  . 

j  =C  Jin ,  (3.90) 


where  Jin is  the  inward  partial  current.  Substituting  this  back  into  equation 
(3.89)  yields: 

-~out  ~zin  -  .  . 

j  =m QiC  Jin+mOEE  ■  (3-91) 


The  outward  partial  current  for  a  given  outward  direction  is: 

Jou,  =  I  "X"  ■ 

ne  out 


(3.92) 


Equation  (3.91)  is  used  to  calculate  the  outward  partial  current: 

J  out  ~  hi  Jin  +  wn  0^OEEhi  • 


neout 


neout 


(3.93) 


The  quantities: 


M0i  -  X  wn  moiC 


neout 


(3.94) 


and 


Je  -  X  Wn{mOEAE )  ) 


neout 


(3.95) 


are  collapsed  coefficients  representing  the  contribution  to  the  outward  partial 
current  in  a  given  direction  from  the  inward  partial  current  and  emissions 
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respectively.  This  process  of  using  the  appropriate  quadrature  weights  to 
integrate  over  a  hemisphere  and  reduce  the  cell  matrices  to  a  single  coefficient  is 
what  I  call  collapsing.  For  equation  (3.95)  the  quantity  JE  represents  the 
outward  partial  current  of  particles  emitted  in  the  cell  that  have  scattered  any 
number  of  times  (0  through  infinity)  before  leaving  the  cell  for  the  first  time.  It  is 
a  known  value  for  the  problem.  For  a  cell  in  slab  geometry,  the  partial  current 
equations  with  the  collapsed  coefficients  are: 

Tl  1  T  Mll  Mlr  T  JL 1  \  JL 

J  Out  _  1V1OI  1V1OI  J  In  ,  JE  /  o 

JR  ~  Mrr  Jr  Jr  ’  1  ' 

J  out  J  L  oi  lvloi  JLT/nJ  \_j e  } 

where  the  superscripts  R  and  L  indicate  the  right  or  left  sides  of  the  cell.  The 
double  superscripts  indicate  the  contribution  to  the  outward  partial  current  from 
the  respective  inward  partial  current,  RL  is  the  contribution  to  the  right  outward 
partial  current  from  the  left  inward  partial  current.  A  similar  relation  can  be 
defined  across  the  spatial  mesh.  Recognizing  that  the  inflow  variables  are 
outflows  of  adjacent  cells,  a  system  of  equations,  Ax  =  b  can  be  set  up  with  the 
emissions  shown  in  equation  (3.95)  as  the  forcing  term  and  the  global  flow 
variables  as  the  unknowns.  Wager  showed  how  the  global  flow  variables 
permutation  resulted  in  a  penta-diagonal  matrix  that  could  be  solved  directly 
(14:  2-47)  using  the  angular  flux  formulation.  The  current  formulation  is  exactly 
analogous.  To  explore  the  efficacy  of  sparse  matrix  methods,  a  Compaq 
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Extended  Math  Library  (CXML)  direct  sparse  matrix  solver  (cxml_dss.f90)  (6: 


11-1)  was  used  to  solve  this  system  of  equations  in  slab  geometry. 


To  distribute  the  partial  current  (or  global  flow  variable)  solution  back  to 
the  detailed  cell  edge  angular  currents,  the  distributions  are  used  as  follows: 


^’L  —  jl  rL 

J  new  Sol'n  ^  "> 

(3.97) 

and 

-R  _  rR  ~fR 

J  new  Sol'n  r 

(3.98) 

A  similar  technique  is  used  in  the  angular  flux  formulation.  This  completes  the 
equations  needed  to  complete  an  iteration  as  described  in  chapter  two. 

C.  Test  Results  In  Slab  Geometry 

1.  Preliminaries:  Measuring  Convergence  Tolerance 
The  symmetric  relative  difference  (SRD)  was  developed  by  Minor  and 
Mathews  (13:  182)  to  determine  when  the  difference  in  the  desired  quantities 
between  iterations  met  the  chosen  convergence  tolerance.  The  relation  for  the 
SRD  is: 


0 


SRD(x,y)  =  l 


2x~y  | 
*\+\y\ 


x  =  y  =  0, 
Else. 


(3.99) 


This  function  returns  a  value  between  zero,  for  values  that  are  exactly  the  same, 
and  two  as  the  limit  for  values  that  are  very  different,  are  of  opposite  signs,  or 
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only  one  of  which  is  0.  Most  often  this  function  is  applied  to  two  arrays  to  find 


the  maximum  SRD  between  corresponding  array  values: 

SRDMax(x,y )  =  Max ( SRD (xt , yt ))  .  (3.100) 

i 

To  check  for  convergence  tolerance,  the  SRD  function  is  applied  to  corresponding 
values  in  two  arrays  (for  two  successive  iterations)  and  the  maximum  value  is 
compared  to  the  convergence  criterion. 

2.  Test  Problems 

The  test  problems  were  a  series  of  single  material  problems  with  a 
reflective  boundary  on  the  left  and  a  vacuum  boundary  on  the  right.  In  the 
series  of  problems,  the  material  was  totally  scattering  with  total  cross  section 
(7  =  1.0  and  an  emission  source  5  =  1.0  uniformly  distributed  throughout  the 
material.  The  cell  size  was  fixed  at  Ax  =  1.0  and  the  number  of  cells  varied  from 
10  to  300  for  the  problems.  The  tests  were  done  using  a  DP8  angular  quadrature 
(16  ordinates).  The  convergence  tolerance  was  10-5  for  the  maximum  symmetric 
relative  difference  (SRD)  in  the  cell  average  scalar  flux  between  two  iterations. 
This  tested  the  performance  of  the  code  for  various  combinations  of  methods  for 
a  series  of  increasingly  larger  problems. 

3.  Edge  Flow  Variable  Formulation 

The  first  series  of  tests  examined  how  DI  performed  with  the  angular  flux 
formulation  as  opposed  to  the  angular  current  formulation  as  discussed  in  chapter 
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two.  The  measure  of  performance  is  the  number  of  iterations  needed  to  reach 


convergence.  The  results  of  these  tests  are  shown  in  table  3.1. 


Table  3.1.  Iterations  to  convergence  for  current  vs.  angular  flux. 


Spatial  Quadrature 

SC 

SC 

LD 

LD 

F  ormulation 

Current 

Flux 

Current 

Flux 

Number  of  cells 

10 

3 

6 

3 

6 

50 

3 

6 

3 

7 

100 

3 

6 

3 

7 

150 

3 

6 

3 

7 

200 

3 

6 

3 

7 

250 

3 

6 

3 

7 

300 

3 

6 

3 

7 

For  these  tests,  discrete  ordinates  sweeping  was  used.  Table  3.1  shows 
that  the  performance  of  the  current  formulation  was  substantially  better  than 
angular  flux  formulation.  The  number  of  iterations  needed  to  reach  the 
convergence  criterion  in  the  current  representation  is  at  most  half  the  number 
needed  for  the  angular  flux  representation.  This  validates  the  choice  of  using  the 
current  formulation  for  additional  implementation  in  XY  -  geometry. 

4.  Coupling  the  Local  Balances 

The  next  tests  examined  the  efficiency  of  different  methods  for  coupling 
the  local  balances  presented  in  chapter  two.  The  test  conditions  were  the  same 
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and  the  tests  were  done  using  the  current  formulation.  The  measure  of 
performance  is  the  total  time  required  to  converge  the  test  problem.  This  was 
measured  using  the  (CPU  time)  intrinsic  FORTRAN  function  and  included  only 
the  computations,  not  file  I/O.  The  results  for  the  step  characteristic  method  are 
presented  in  figure  3.1. 


Figure  3.1.  The  problem  solution  time  versus  number  of  cells  for 
different  cell  flow  coupling  methods  with  step  characteristic. 

Figure  3.1  shows  that  the  discrete  ordinates  sweep  is  the  most  efficient 
method  for  coupling  the  local  balance  among  cells  for  this  single  processor 
implementation.  The  Red/Black  method  performance  was  less  efficient  than  the 
discrete  ordinates  sweep  but  still  an  improvement  over  the  local  balance  sweeping 
method. 
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The  results  of  the  same  test  for  linear  discontinuous  with  the  current 


representation  are  shown  in  figure  3.2. 


Number  of  Cells 

Figure  3.2.  The  problem  solution  time  versus  number  of  cells  for  different 
cell  flow  coupling  methods  with  linear  discontinuous. 


The  results  for  linear  discontinuous  method  confirm  the  results  shown  by  the  step 
characteristic  method.  A  similar  timing  test  was  run  for  a  two  material  problem, 

with  an  absorbing  ( c  =  —  =  0.3  )  region  and  a  scattering  ( c  =  0.3  )  region.  The 

total  times  were  again  consistent  with  the  previous  series  of  one  material 
problems:  the  discrete  ordinates  sweep  method  had  the  fastest  time  and  local 
balance  sweeping  the  slowest.  The  overall  efficiency  of  the  discrete  ordinates 
sweep  method  makes  it  a  good  choice  for  implementation  in  XY  -  geometry 
because  it  will  be  demonstrated  on  a  serial  machine.  For  problems  large  enough 
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to  need  parallel  implementation,  the  Red/Black  method  may  be  the  method  of 
choice. 
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IV.  Implementation  in  XY  -  Geometry 


This  chapter  presents  the  derivations  for  the  local  balance  problems  for  both 
zeroth  and  first  spatial  moments  problems  using  the  current  representation.  The 
global  balance  problem,  subsequently  called  the  partial  current  problem,  is  also 
presented. 

A.  Zeroth  Spatial  Moment  Methods  Distribution  Iteration  Derivation 

The  zeroth  spatial  moment  methods  are  an  extension  of  the  method 
presented  in  slab  geometry  in  chapter  three.  The  general  representation  used  in 
chapter  three  is  changed  to  explicitly  account  for  the  contributions  for  each 
cardinal  direction,  even  though  in  some  cases  the  contribution  is  zero. 

The  desired  form  is  to  assemble  the  equations  for  zeroth  spatial  moment 
methods  in  a  relation  that  gives  the  cell  face  outgoing  currents  in  terms  of  the 
cell  face  inward  currents  and  cell  emissions.  A  general  method  for  all  zeroth 
spatial  moment  methods,  such  as  step  characteristic  (SC)  or  weighted  diamond 
difference  (WDD),  is  presented  here. 

The  form  of  the  system  of  equations  is: 


- *•  out  —  in  —  — 

j  -  K0/  j  +  K osS  +  K0eE  > 

(4.1) 

—in  —  — 

(4.2) 

and 

S  =  Is  V  ■ 

(4.3) 
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In  equation  (4.3),  the  scattering  matrix  Z5  has  the  same  representation  as 


in  equation  (3.3).  The  vectors  are  defined: 


-out 

j 


-out 

J L 


—out 

Jr 


—out  —out 

Jt  J  B  > 


(4.4) 


— in 

j 


—in  —in  —in 

j l  Jr  Jt 


(4.5) 


and 


S  =  SA,  (4.6) 

E  =  E4.  (4.7) 


In  equations  (4.4)  and  (4.5)  the  directions  for  the  sub-vectors  are  given  by  the 


capital  subscript,  for  example  L  for  left.  The  matrices  are  defined: 


Ko/  - 


1 \{LL 

*^OI 

ts  LR 
*^OI 

1 {LT 
K OI 

K 

1 rRL 

*^OI 

URR 

*^OI 

ICRT 

*^OI 

k; 

\zTL 

*^OI 

ktr 

ktt 

K 

is  BL 

Ko/ 

kbr 

K OI 

kbt 

*^OI 

K 

LB 

01 

RB 

OI 

TB 

OI 

BB 
OI  . 


Kos  - 


L 

OSA 

R 

OSA 

T 

OSA 

B 

OSA . 


K 


OE 


K 

K 

K 

K 


L 

OEA 

R 

OEA 

T 

OEA 

B 

OEA . 


K  AI  = 


K 


AI 


K 


R 

AI 


K 


AI 


K 


AI 


K 


t//S 


K 


ASA  1 


(4.8) 


(4.9) 


(4.10) 


(4.11) 

(4.12) 
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and 


(4.13) 


-  K aea  ■ 

In  equations  (4.8)  through  (4.11)  the  sub-matrices  are  diagonal  matrices  similar  to 
those  used  in  slab  geometry.  In  equation  (4.8)  the  sub-matrices  are  matrices  that 
give  outgoing  currents  from  inward  currents,  hence  the  subscript  01,  and  the  two 
directions  in  the  superscript  correspond  to  the  outward  direction  from  the  inward 
direction.  For  example,  K0/  is  the  sub-matrix  that  gives  the  right  outward 
currents  from  the  left  inward  currents.  In  equations  (4.9)  through  (4.10),  the 
superscript  directions  correspond  to  the  outgoing  direction  and  the  subscripts 
have  the  same  meaning  as  in  slab  geometry.  For  example,  K OSA  is  the  matrix 
giving  the  current  out  the  top  from  the  average  scatter,  K qEA  is  the  matrix 
giving  the  current  out  the  left  from  the  average  emissions,  and  K AI  is  the  matrix 
giving  the  average  flux  from  the  right  inward  current. 

Similar  to  the  slab  geometry  case,  scatter  is  eliminated  from  equations  (4.1)  and 
(4.2): 

V  =  (I  -  Kw  Zs)-'[Kr,>'"  +  Kr«£] .  (4.14) 


Equation  (4.14)  can  now  be  used  to  eliminate  scatter  from  equation  (3.50).  The 
resulting  equation  is: 


]-• = <K„ + kos  zs(i  -  Krs  zsr‘  k  rI)r + 

(K0£  +Kos  Es(I  -  Krs  Z  s)-‘  K„)£. 


(4.15) 


This  equation  can  be  used  to  calculate  the  cell  outgoing  currents  from  the  cell 
inward  currents  and  cell  emissions.  Looking  at  the  terms  in  each  parentheses, 
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the  first  matrix  is  the  uncollided  contribution  to  the  respective  outgoing  current 


from  the  respective  incoming  current  as  modeled  by  the  spatial  quadrature 
method.  The  product  or  second  term  represents  the  respective  incoming  current 
that  contributes  to  the  respective  outgoing  current  after  all  scattering  takes 
place,  again  as  modeled  by  the  spatial  quadrature  method. 

In  addition,  the  final  matrix  represented  by  the  sum  in  each  parentheses, 
only  needs  to  be  calculated  once  for  each  combination  of  cell  size  and  material. 


The  final  matrices  for  equation  (4.15)  can  be  expressed  as  a  matrix  equation: 

- ~out  -in  —A  . 

j  =m oij  +m qee  >  (4-16) 

where  mOI  is  a  matrix  that  gives  outward  currents  from  inward  currents  and  mOE 
is  a  matrix  that  gives  outward  currents  from  emissions.  The  current  vectors  and 


emission  vector  with  the  respective  sub-matrices  for  each  direction  is: 


J  Out 

J  Out 

Tout 

i 

b-i 
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OI 


m 

m 

m 

m 


01 

LR 

Ol 

TR 
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Ol 
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OI 

BB 

OI 


Tin 

KJA 

Ifn 

+ 

<E JA 

Tin 

<eaEA 

Ji\ 

<eaEA_ 

(4.17) 


Here,  the  sub-matrices  represent  the  outward  contribution  from  the  inward 


current  after  any  number  of  scatters.  The  emissions  vector  represents  the  forcing 
term  for  this  system  of  equations.  Later,  the  relation  between  the  inward  and 
outward  currents  will  be  used  to  set  up  the  partial  current  problem  across  the 


spatial  mesh.  The  next  section  will  show  how  to  calculate  the  values  for  the 
diagonal  sub-matrices  in  equations  (4.8)  through  (4.13). 
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1.  Step  Characteristic  Transport  Coefficients 


The  first  spatial  method  implemented  in  XY  -  geometry  was  step 
characteristic  (SC).  This  method  was  chosen  for  its  relative  simplicity,  and  as  a 
way  to  validate  the  extension  to  XY  -  geometry  before  attempting  other  more 
complicated  spatial  methods.  Miller  (10:  21)  presents  the  cell  equations  in  the 
angular  flux  representation  using  the  exponential  moment  functions.  The 
derivation  for  the  current  equations  is  in  appendix  B. 

Top 


Bottom 

Figure  4.1.  Rectangular  cell  for  implementation  of  the  zeroth  spatial 
moment  methods. 


For  the  cell  in  figure  4.1  showing  ordinate  n  out  the  top  face,  JUn  and  l]n  are  the 
direction  cosines  along  the  x  and  y  axis  respectively  from  the  angular  quadrature 
set,  the  optical  thickness  along  an  ordinate  in  the  y  and  x  direction  is: 

crAy 


"yn 


fcl  ’ 


(4.18) 
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(4.19) 


g.yH  _|/6,|4 y 

£xn  ' 


(4.20) 


These  equations  are  used  for  the  spatial  quadratures  in  XY  -  geometry.  For  a 
rectangular  cell  as  shown  in  figure  4.1,  the  equations  for  the  outgoing  currents  in 
ordinate  n,  j*p  and  jn"ght ,  in  terms  of  the  incoming  currents  j„bot,om  and  jr‘efi  , 
scattering  within  the  cell,  SAn ,  and  emissions,  EAn ,  for  ordinate  n  with 
rjn  >  /Jn  >  0  and  an  <  1  are: 

=n«£,.)^ 

\Vn\  (4.21) 

+  anMx  (£y,n  )]S?  +  Av[(l  -  an  )M0  (eyn )  +  anMx  )]E^ , 

j  =  (e  +hlYMl(£,,  )SA  +  1 )EA ,  (4.22) 

I  j-j  y  II  y  11  y  n 

\'/n\  \'ln\  'In 


Yn  =T—^nMo(£yn)jnLeJt+r-][(l-^n)Mo(£yn)  +  (XnMl(£yn)\jnBOtt0m  + 

|  Mn\  \Vn\ 

— [(1  -  an  )Ml  (£yn )  +  anM2  )]S?  + — [(1  -a„  )M1  (£},n )  +  anM2  )]Ef , 

'In  'ni 


(4.23) 


where  M0(£v  ),  M{(£v  )  and  M2(£v  )  are  the  exponential  moment  functions 


defined  in  equation  (3.9).  Because  the  ratio  of  the  direction  cosines  is  frequently 
used,  let: 


(4.24) 


65 


The  transport  coefficients  are  the  values  used  to  build  the  diagonal 
coefficient  matrices  described  previously.  Equations  (4.21),  (4.22)  and  (4.23)  can 
be  written  as: 


Top  _  TB\  ■  bottom  ,  /  r  TL\  ■  left  ,  si  TR\ 

Jout,,  ~  \K-OI  )nJin„  '  \^OI  hi  Jin„  '\^OlhiJin. 


right 


+ 


m  )»  im, ,0P  +  (toSA  >„  V  +  d&EA  )„  E„  , 

right  /  /  RB  x  •  bottom  ,  /  j  RL  \  ■  left  ,  /  r  RR  \  ■  right  , 
Jout„  =  (kOl  )nJin„  +(kOl)nJinn  +(kOI>nJinn  + 

Ki  )J,n.,op  HkosA)AA  +(k%EA)nE„A, 


(4.25) 


(4.26) 


and 


(4.27 

(b/  )„  J„fP  +  (K<SJ  )„  Sf  +  (kAEA  )„  E,f 
This  form  is  a  variation  of  equation  (3.50)  giving  the  outgoing  currents 
from  the  inward  currents,  scattering  and  emissions  for  an  ordinate.  As  was  done 


in  slab  geometry,  the  k  vectors,  which  are  used  to  form  the  diagonal  matrices 
K  =  D(fc)  ,  can  also  be  found  by  inspection  of  equations  (4.21),  (4.22)  and  (4.23) 
Thus  the  k  vectors  for  the  top  outward  angular  currents  are: 


(kSj)n=(\-an)e-£>»  , 

(4.28) 

(kOl)n=  *nanM0(£y„)i 

(4.29) 

ikToRI)n=(kTJI)n=  0, 

(4.30) 

and 

&OSA)n  =  ( kOEA)n  =  A.v[(l  -  OSn)M0(£yn  )  +  (Xn M ^  )]  . 

(4.31) 

The  k  vectors  for  the  right  outward  angular  currents  are: 


66 


(uRB\  -  M^£yJ 

\KOI  )n  -  > 

^ n 

(4.32) 

(k Ol  )n  =  (kof  )n  =  (koi  )n  =  0  , 

(4.33) 

and  (kosA )n  =  (k(jFA )n  =  — —  M\(£yn )  ■ 

~n 

(4.34) 

The  k  vectors  for  the  average  angular  flux  are: 

u  1 

(kAI)n  =  ^[(1  -  an)M0(£yn )  +  anMx{£yn )] , 

(4.35) 

(kLAl)n  =  1  ^nM0(£yn), 

\  Mn  \ 

(4.36) 

(4>„ =(b/)„= o. 

(4.37) 

Av 

and  (kASA  )n  =  (kAEA  )n  =  •  [(1  an  )Ml  (£yn )  +  anM2  (^ )] . 

'In 

(4.38) 

Each  ordinate  is  evaluated  to  determine  the  outgoing  face  and  the 
respective  transport  coefficient.  Not  shown  in  figure  4.1  are  ordinates  exiting  the 
right,  bottom,  or  left  cell  edges  instead  of  the  top  edge;  however,  the  same  basic 
relations  are  used.  For  these  cases,  an  x-y  reversal,  a  right-left  exchange  or  a 
top-bottom  exchange  are  used  where  appropriate.  In  all,  there  are  30  transport 
coefficients  to  find  for  each  ordinate  of  which  only  11  are  nonzero. 

These  transport  coefficients  are  used  to  build  the  diagonal  matrices  listed 
in  equations  (4.8)  through  (4.13).  Once  the  diagonal  matrices  are  calculated,  the 
final  matrices  described  in  equations  (4.15)  and  (4.17)  can  be  constructed.  For 
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the  DI  method,  these  matrices  will  be  calculated  one  time  for  each  material 
(assuming  the  spatial  mesh  is  uniform). 

2.  Weighted  Diamond  Difference  Transport  Coefficients 

Another  zeroth  spatial  moment  method  that  was  implemented  was 
weighted  diamond  difference  (WDD).  This  positive  method  was  chosen  to 
demonstrate  the  ease  of  adding  other  zeroth  spatial  moment  methods  and  to 
compare  with  published  results.  Azmy  (3:  215-216)  presents  the  angular  flux 
formulation  of  the  WDD  method,  which  is  changed  to  the  current  formulation  in 
appendix  B.  For  the  same  rectangular  cell  shown  in  figure  4.1,  the  equations  for 
the  outgoing  currents  jj0’’  and  jnnght ,  in  terms  of  the  incoming  currents  j nbollom 
and  jnlefl  scattering  within  the  cell  S)'1  and  emissions  A4 ,  for  ordinate  n  with 
r/n  >  jUn  >  0  and  an  <  1  are  presented  in  Appendix  B.  To  avoid  bad  numerical 
conditioning,  the  WDD  equations  are  cast  in  terms  of 
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(4.42) 


K  = 


SX"  £ 

uoutc'x, 


"  c 

uoutcyn 


8X"  £ 

. OUtc'X . 


+  8^n  £  +  8X"  £ 

'  uoutt'y  '  uoutt'x, 


%y n  c 
uoutc’yn 


Determining  the  values  of  the  k  vectors  used  to  form  the  diagonal  matrices 
K  =  D(£)  is  done  the  same  way  that  was  used  for  the  SC  quadrature.  The  k 


vectors  for  the  top  outward  angular  currents  are: 


(1TB,  _  K  8\n 

\kOI  )n  eV  7  eV  ’ 

^Out)  £yn  5Out 

(4.43) 

(]'TL\  _  hnK\ 

k  r 

(4.44) 

(4.45) 

and 

nT  \  nT  \  KYln\ 

\kOSA)?i  ~\kOEA)?i  ~  _v 

(4.46) 

The  k 

vectors  for  the  right  outward  angular  currents  are: 

l£RB\  -  K\Rn\ 

t  01  )n  §yn  I  I  > 

^ Out  ^ Out  ^ yn  Y/\ 

(4.47) 

/  j_RL\  hn  ^ In 

(4.48) 

(O„  =  (O„  =  0, 

(4.49) 

and 

(]R  \  -(/CR  \ 

y^OSASn  y^OEAJn  ?x 

°Out® 

(4.50) 

The  k 

vectors  for  the  average  angular  flux  are: 

ii 

R 

(4.51) 
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(A)„  = h- ,  . 

^ Out  \Mn  \  ^ xn 

(4.52) 

(^Al)n  =(^Al)n  =0  > 

(4.53) 

and 

h n 

(^ASa)h  = (k AEA)n  =  ~  ■ 

(4.54) 

Unlike  characteristic  methods,  the  WDD  equations  need  not  treat 
Ax  /  \n\  <  Ay  /  |?y|  differently  than  Ax  /  |/x|  >  Ay  /  \r]\ .  For  convenience  in  sharing 
code,  I  use  the  equations  for  the  case  in  figure  4.1  to  fill  the  WDD  matrices  in  the 
same  way  as  I  described  above  for  SC.  Again,  for  the  DI  method,  these  WDD 
matrices  (or  any  other  linear,  zeroth  spatial  methods)  use  the  same  solver 
algorithm  as  SC. 

B.  Derivation  of  First  Spatial  Moment  Methods  Distribution  Iteration 

Similar  to  the  zeroth  spatial  moments  derivation,  it  is  desirable  to 
assemble  the  discrete  ordinates  system  of  equations  in  a  form  that  gives  the  cell 
outgoing  currents  in  terms  of  the  cell  inward  currents  and  cell  emissions. 

However,  unlike  the  zeroth  spatial  moment  methods,  there  is  the  addition  of  the 
first  spatial  moment  of  the  current  along  edges  0  and  first  spatial  moment  of  the 
scattering  sources  in  each  dimension  to  consider.  The  equations  could  again  be 
cast  into  the  general  form: 


—  out  -~in  -*  - 

j  =  K-oi  j  +  +  K oeE  , 

(4.55) 

_  —in  —  — 

¥  =  k^/7  +  +  ^yeE  > 

(4.56) 

and 

S  =  Ts¥, 

(4.57) 
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where  the  angular  current  vectors  include  the  first  spatial  moment  of  the  current: 


and 


—in 

j 


(4.58) 


-out 

j 


-’OUt  -’OUt 

j  0 


(4.59) 


This  approach  was  not  used;  instead,  the  spatial  moments  of  the  angular  currents 
are  explicit  in  the  system  of  equations.  This  simplifies  the  indexing  for  the  code 
and  allows  use  of  the  same  routine  for  the  partial  current  problem  as  for  the 
zeroth  spatial  moment  quadratures.  This  routine  is  presented  later  in  this 
chapter.  The  cell  system  of  equations  for  the  first  spatial  moment  spatial 
quadratures  in  general  can  be  written: 


— out  —in  —in  —  — 

j  =  KOlJ  +KO0#  +KOS‘S  +  K0£'jE'’ 
eout  =  K  eij  +  K00oin + K0Ss  +  K  0Ee, 

—  —  in  —in  —  — 

¥  =  ^Vl]  +Kp6>^  +KV/SS  +  KiyEE> 


and 


s  =  'Ly/- 


The  vectors  for  equations  (2.2)  through  (2.4)  are  defined  as: 


- out 


—out  —out  —out  —out 
Jl  Jl  j  L  j  L 


.  nT 

—in  —in  —in  —in 

Jl  Jl  Jl  J  l 


-out 

0  = 


—out  -out  -out  -out 

eL  Ol  9l  9l 


~\T 


e 


—  in  —in  —in  —in 

9l  Ol  0l  Ol 


(4.60) 

(4.61) 

(4.62) 

(4.63) 

(4.64) 

(4.65) 

(4.66) 

(4.67) 
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(4.68) 


and 


(4.69) 

(4.70) 


The  matrices  for  equations  (2.2)  through  (2.4)  follow  the  same 


methodology  as  slab  geometry.  The  matrices  for  equation  (2.2)  are  defined: 


and 


K07  - 


K  o*  = 


1 \{LL 

Ko/ 

1 7-LR 
^07 

^07 

\iRB 

^07 

1 {RL 
Ko/ 

URR 

krt 

K07 

1^75 

K07 

1 {TL 
k07 

ktr 

K07 

ktt 

K07 

^75 

K07 

\zBL 

Ko/ 

kbr 

K07 

kbt 

K07 

K07 

1 {LL 
KO0 

is  LR 
KO0 

KO0 

\zLS 

K06> 

1 -s-RL 
KO0 

krr 

*^oe 

^77 

K06> 

^75 

K0<9 

1 {TL 

K oe 

ktr 
K oe 

tt-TT 

*^oe 

^75 
K 0(9 

\zBL 

K oe 

kbr 

Ko<? 

kbt 

KOi9 

^55 

KO0 

*VaS>1 

^OSX 

SO 

* 

1^7 

KOi>l 

\iR 
l'-  osx 

i^7 

^057 

^ OS  A 

^ OSX 

v'T’ 

^057 

\^B 

.0571 

mb 

^ osx 

1 sB 

^057  J 

^  07,1 

^OEX 

K0£Y 

1^7 

^07/1 

KR 

OF.X 

1^7 

^ OEY 

,v07/( 

i yB 

OF.X 

^ 077 

lv07/l 

MB 

IV  OF.X 

1^7 

^077. 

(4.71) 


(4.72) 


(4.73) 


(4.74) 


The  notation  is  similar  to  the  zeroth  spatial  moment  method  notation.  For 
example,  the  diagonal  sub-matrix  represents  a  matrix  that  gives  bottom 
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outward  currents  from  left  inward  currents.  The  diagonal  sub-matrix  K  'osx  gives 
the  left  outward  current  from  the  x-moment  of  scatter  in  the  cell  S  .  The 


diagonal  sub-matrix  K ROEA  gives  the  right  outward  current  vector  from  the 

—  A 

average  emissions  in  the  cell  E  .  The  higher  moments  for  the  emissions  vector, 
as  for  the  first  moments  of  emissions  in  slab  geometry,  represent  the  higher 
moments  of  scatter  from  other  energy  groups  in  a  general  representation.  For  the 
mono-energetic  problems  used  in  this  research,  these  vectors  and  matrices  are  not 


used.  The  matrices  for  equation  (4.61)  are  defined: 


K  01  = 


K  ee  = 


1 ifLL 

*^0I 

1 {LR 

*^0I 

\zLT 

K6>/ 

is  LB 
K<97 

1 rRL 
*^0I 

! v-RR 

*^0I 

\CRT 

*^0I 

krb 

*^0I 

1 {TL 

*^0I 

KTR 

K<97 

KTT 
K 01 

ktb 

*^0I 

\zBL 
K 01 

kbr 

*^0I 

kbt 

Ki97 

KBB 

*^0I 

1 fLL 
K 00 

1 sLR 
K 00 

KLT 

\iBB 

K<9<9 

URR 
K 00 

tzRT 

^00 

URB 

K<9<9 

1 {TL 
K 00 

ktr 
K 00 

ktt 

K00 

ktb 
K 00 

1 fBL 
W00 

kbr 
K 00 

1 zBT 
K 00 

Kbb 

K00 

K 


os 


K 

K 

K 

K 


L 

0SA 

K  OSX 

K 

R 

0SA 

K  OSX 

K 

T 

0SA 

1 (T 

K  OSX 

K 

B 

0SA 

IV  OSX 

K 

L 

osy 

R 

osy 

T 

osy 

B 

0SY . 


and 


(4.75) 


(4.76) 


(4.77) 
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KffEA 

^0F,X 

K GEY 

KffEA 

^GEA 

ts  T 
W0F,X 

T 

^ 0EY 

^GEA 

1 ~sB 
^GEX 

| irB 

^ GEY 

(4.78) 


Again,  the  notation  is  similar  to  the  previous  notation  used  for  equation  (2.2). 
Here,  the  outgoing  current  symbol  O  is  replaced  with  the  outgoing  edge  current 
moment,  9  so  the  diagonal  sub- matrix  represents  a  matrix  that  gives 

bottom  outward  edge  current  moment  from  left  inward  currents.  The  matrices  for 
equation  (2.3)  are  defined: 


(4.79) 


1 

* 

’ — 1 

kr 

K/l/ 

K  TAi 

1 

* 

K -x, 

*^XI 

K  Txi 

Kv/ 

1 

* 

K  YI 

Kyi 

1 sB 

K  YI 

K  LAe 

kAG 

K  TAg 

K X6 

t^R 

K XG 

KTxe 

K YG 

j^R 

K  tY9 

(4.80) 


KASa  K asx  K asy 
=  K XSA  K xsx  KXSY  j 


(4.81) 


*^AEA  ^ AEX  *^AEY 

Ky/E  =  KXEA  KXex  KXey  ’ 


(4.82) 


I,  0  0 

0  I,  0  . 
0  0  X, 


(4.83) 
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Similarly,  the  diagonal  sub-matrix  K74/  gives  the  contribution  the  average 

—  .4 

angular  flux  within  a  cell  !//  from  the  top  inward  current  vector,  the  diagonal 

sub-matrix  gives  the  contribution  the  x-moment  angular  flux  within  a  cell 
— x 

iff  from  the  right  inward  edge  current  moment  vector,  the  diagonal  sub-matrix 
K  XSY  gives  the  contribution  the  x-moment  angular  flux  within  a  cell  !//  from  the 
y-moment  of  scatter  vector,  and  the  sub-matrix  X<,  is  the  scattering  matrix 
defined  in  equation  (3.5). 

Equation  (4.63)  is  substituted  into  equation  (4.62)  to  eliminate  scatter. 

This  gives: 

^  =  (I-KrSE)-1A:f,;y'+(I-KrSZ)-|K(M,r'+(I-Kf/SS)-1KrJ;I.  (4.84) 
Equation  (4.84)  is  then  substituted  into  equations  (4.60)  and  (4.61)  to  again 
eliminate  scatter.  The  final  equations  are: 

7"  =  (Kc;  +  Kos  HI  -  K^s  S)-1  K r,)7  +(Kot)+  Kos  HI  -  E)"1  K ^0"  + 

(K0£  +  Kos  2(1  • -  K„s  »_1Kr£)I, 

and 

7"  =(Kw+KMHI-K^Sr1K^)7"  +(K«,  +  K9sZ(I-KrSZr1K^)§‘"  +(J 
CK»£  +  Kffi  HI  -  K„s 

As  was  done  for  slab  geometry,  the  code  implementation  used  an  equivalent 
elimination  of  the  scattering  sources,  in  which  a  block  forward  elimination  and 
back  substitution  produces  the  inverse  matrix  in  equations  (4.84)  through  (4.86), 
which  is  shown  in  appendix  C. 
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As  with  the  zeroth  spatial  moment  methods,  the  final  matrix  represented 
by  the  sum  in  each  parentheses  for  equations  (4.85)  and  (4.86),  only  need  to  the 
calculated  once  per  material  and  cell  size.  The  matrices  for  these  equations  can 
be  expressed  as: 


-out 


*1 out 


-tin  —  A 

,0  +m oeE  , 

(4.87) 

-tin  —  A 

e  +m  eEE  . 

(4.88) 

and 

Here  m()l  is  a  matrix  that  gives  cell  outward  currents  from  inward  currents, 
m0£is  a  matrix  that  gives  cell  outward  currents  from  emissions,  mw  is  a  matrix 
that  gives  cell  outward  current  edge  moments  from  inward  currents,  and  m gg  is  a 
matrix  that  gives  cell  outward  current  edge  moments  from  inward  current  edge 
moments.  The  current  vectors,  edge  moment  vectors  and  emission  vectors  with 
the  respective  sub-matrices  for  each  direction  is: 


Ti 

J  Out 

l 

3 

8  S 

m" 
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m 

R 

J  Out 

< 

< 
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m 

J  Out 
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m 

T  B 
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m" 

nC 

m 

m 

m 
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LL 

oe 

RL 

oe 

TL 

oe 

BL 

oe 


m 

m 

m 

m 


LR 

oe 

RR 

oe 

TR 

oe 

BR 

oe 


m 

m 

m 

m 


LT 

oe 

RT 

oe 

TT 

oe 

BT 

oe 


m 

m 

m 

m 


TL  " 

J  In 

T  R 

J  In 

+ 

TJ 

J  In 

_ 

_J  In  _ 

~  — L  “ 

din 

L  t?A 
111 OEA1^ 

—  R 
din 

+ 

<eaEA 

—  T 

din 

<eaEA 

- 

-* B 
din 

<eaEA_ 

(4.89) 


A  similar  system  can  be  constructed  for  the  outward  current  edge  moments  for  a 
cell: 
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1 

0 

K 

_ 1 

r 

0  Out 

-T 
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-~B 

0  Out 
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_ LT 
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m9EA E 
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(4.90) 


Again,  the  sub-matrices  represent  the  outward  contribution  for  the  respective 
vectors  from  the  inward  vector  after  completing  any  number  of  scatters.  The 
emissions  vector  represents  the  forcing  term  for  these  systems  of  equations. 

These  equations  will  be  used  to  set  up  the  partial  current  problem  across  the 
spatial  mesh.  The  next  section  will  show  how  to  calculate  the  values  for  the 
diagonal  sub-matrices  in  equations  (4.71)  through  (4.82). 

1.  Linear  Characteristic  Transport  Coefficients 

The  first  method  implemented  was  linear  characteristic  (LC)  which  was 
initially  developed  by  Alcouffe  et  al.  in  1979  (11:  24).  This  first  spatial  moment 
method  was  chosen  as  an  extension  of  SC  and  to  show  the  implementation  of 
first  order  spatial  methods.  Miller  (12:  23)  also  provided  the  LC  cell  equations  in 
the  angular  flux  representation  using  the  exponential  moment  functions.  The 
derivation  for  the  cell  current  equations  are  presented  in  Appendix  D.  The 
process  of  determining  the  equations  for  the  k  vectors  used  to  form  the  diagonal 
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matrices  K  =  D(£)  is  analogous  to  the  procedure  used  for  SC.  The  k  vectors  for 


the  top  outward  angular  currents  are: 

(k™)n  =  (l-an)e-£r»  ,  (4.91) 

(ka)n=TnanM0(£yn),  (4.92) 

(4=^0-^,  (4.93) 

(koe)n  =Tna„[2M1(£yn)-M0(£yn)] ,  (4.94) 


(kosA)n=(koEA)n=anM-(l-an)M0(£yn)  +  (\-2an)Ml(£yn)  +  anM2(£yn)\,  (4.95) 

(kosx)n  =^y\-^-(Xn)MQ(£yn)  +  {\-2an)Ml(£yn)  +  anM2(£yn)\ ,  (4.96) 

and  (koSY)n=Av[-(l-an)M0(£yn)  +  (2-3an)Ml(£yi)  +  anM2(£yi)\.  (4.97) 

The  k  vectors  for  the  right  outward  angular  currents  are: 


Tn 

(4.98) 

(koe  )n  =  — [(1  -  2  an  )M0  (eyn )  +  2anMx  )] , 

(4.99) 

( koi  )„  =  (kf^  )„  =  0  , 

(4.100) 

(kosy)n  =  ~\_k4  2  (£yn  )~MX  (£y,n )] , 

(4.101) 

(kosx)n  =  —[0  - 2«;I)A7, (£y„ )  +  2anM2(£yn )] , 

Tn 

(4.102) 

( koSA)n  =  (k()EA)n  =  — Mx(£y  )  . 

(4.103) 

The  k  vectors  for  the  top  outward  first  moment  of  the  angular  currents  are: 

(C)„=3«„(l (4.104) 
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(kT0i  )n  =  3Tnan[(2an-\)M0(eyn)-2anMx(eyn)] , 


(4.105) 


(kT0B0)n=(\-3an+2an2)e~£^ 


(kTee)n  =2Tnan[(l-2an)M0(£ y)  +  (6a„ -2)Ml(s  )-4anM2(£  )] , 


(keSA  )n  =  (keEA  )n  =  3anAy[(l  -  an  )M0  {eVn  )  +  (2 an- 1  )Ml  (eK )  -  ocnM2  (£yy )] , 


ov 


ov 


(klsx  )n  =  Ay[(l  -  3 a„  +  2a,3  )M0  (£y,n )  +  (3a n  -  6a3)Mx  )  + 

6  a3M2  (£yn )  -  2a,3  M3  (£}>n )] , 


and 


(klsY  )n  =  3«;IAv[-(l  -  an  )M0  (£,„ )  +  (3  -  4a, ,  )Ml  )  - 

(2  -  5an  )M2  (£yn )  -  2  anM3  )] . 


The  k  vectors  for  the  right  outward  first  moment  of  the  angular  currents 

(kSf)n=—[M0(ey)-2M1(ey)], 


( kw)„  =  —[(1  - 2 an)M0(£  )  +  (6 an  - 2 )Ml(e  ) - 4 anM2(£  )] , 


(*«%=(*#)„  =0, 


C k$SA)n=(kgEA)n  =—[M1(£  )-M2(£  )\  , 


( kesxX  =  —  [(1- 2an)Mx  (£yn)-<}- 4an)M2(£yn)- 2anM3(£yn )] , 

Tn 


and 


(kjjsY)n  -  —\~^M\(£yn  )  +  k>M2(£yn  )  -  2  M3(£yn  )] . 
Tn 


The  k  vectors  for  the  average  angular  fluxes  are: 


D  1 

(kA/ )„  =  [(1  -  a„  )Mq  (£yn )  +  anMx  )] , 


(4.106) 

(4.107) 

(4.108) 

(4.109) 

(4.110) 

are: 

(4.111) 

(4.112) 

(4.113) 

(4.114) 

(4.115) 

(4.116) 

(4.117) 
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K  )„=prMi  (£„),  (4.118) 

|  An  | 

(kBAo)n  =  -p\[-(l - an )M0 (e  )  +  (1  - 2 an )Ml (e  )  +  anM2 (e  )] ,  (4.119) 

yin  | 

(kAe)„  =  pr  (£,_)],  (4.120) 

W.  =  (kAEA)„  =  )‘  [H-ffjMilf,.,  )  +  (4.121) 


( W )„  =  %T[-d  -  ■ (  V )  +  ■ (!  ■ -  2«„  )M,  (£r> ) + a„M}  (e^ )] ,  (4.122) 

yin  | 

and  (W)„  =riHl-«nW1(£ k)  +  (1-2«,;W2(^ y)  +  OCnM3{£ )] .  (4.123) 

'yy  x  ll  -'ll  s  tl 

yin  | 

The  k  vectors  for  the  x-moment  of  the  angular  fluxes  are: 


(*£)„  =^[(1  — «/;)M0(^)  +  (2a,J -\)Ml(e.,)~a„M2(ey ,  )] , 


(*&  )n  =  ■ ^[(2^  -  1)M!  )  -  2 anM2  )] , 

An 


(4.124) 

(4.125) 


(*!*)«  =  ^[(l-3«n  +  2a„3)M0(£- „)  +  (3orw -6a3)Mx{£  )  + 


I  An  | 

6a3M2  (£  )  -  2a3M3  (e  )], 


(4.126) 


(*!»)»  =f)[(l- 2ff„  )  -  (4«„  -  mi(e>A )  -  2anM,{eyn )] , 

|An  | 

(*!/)„  =  “  “»>M o(ey. ) +  (2«„  -  DM,  (£ )  -  a„M2 (e  )] , 

y/n  I 

(kXSA  )n  =  (kXEA  )n  =  IX1  ~  <*n  )M\  (£y„  )  +  (2«n  “  DM2  )  -  a„M3  )]  , 

‘In 


(4.127) 

(4.128) 

(4.129) 
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(4.130) 


and 


(kxsx)n  =  “  3a»  +  2an  )M\  (£yn  )  +  ( 3an  ~  )M2  (£yn  )  + 

yin  | 

6a„3M3  (£yn )  -  2  a„3M4  (e^ )], 


ikxsr  )n  =  -  «„  )M,  (£„ )  +  (2  -  3a, ,  )M2  )  - 

(1  -  3or„  )M3  )  -  a„M4  (eyn )] . 


(4.131) 


The  k  vectors  for  the  y-moment  of  the  angular  fluxes  are: 

(ky,)n  =^-\(l-an)M0(£yn)  +  (3an-2)Ml(eyn)-2anM2(£yn)\ , 

'In 


(kyi)r 


3 On 

Vn 


\M\(£yn)-M2(£vJ\, 


«™)»=?n-Hl-«»)*'oVr,)  +  (3-4a„)M1(£  )+ 

I'M 

(5a„  - 2)M2(£  ) - 2 anM3(£  )], 


(4.132) 

(4.133) 


(4.134) 


(*y*)«  =  ^[-3M0  )  +  6Mj  (£yn  )  -  2M2  (£yn )], 

(kySA)n  =  (kYEA  )„  =  I [fl "  Wl  (^„  )  +  (2^  "  [)M2  (£yn  )  "  ^nM3  (£y„  )]  , 

few )„  =  ^[-(1  -  a, )M,  (e  ) + (2  -  3 a„)M2 (e  )  - 

*  H  *  n 

y  In  | 

(1  -  3or„  )M3  (£yn )  -  anM4  )] , 

(kysy)n  =  r7[-3(! - an)Mx(£  )  +  (6-9 an)M2(e  ) - 
and  yin  \ 

(2  -  8  an  )M3  (£yn )  -  2  anM  4  )] . 


(4.135) 

(4.136) 


(4.137) 


(4.138) 


As  with  the  zeroth  spatial  moment  methods,  each  ordinate  is  evaluated  to 
determine  the  outgoing  face  and  the  respective  transport  coefficient.  Ordinates 
exiting  the  right,  bottom  or  left  cell  edges  instead  of  the  top  edge  use  the  same 
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basic  equations  with  an  x-y  reversal,  a  right-left  exchange  or  a  top-bottom 
exchange  where  appropriate. 

2.  Linear  Discontinuous  Transport  Coefficients 

The  next  spatial  method  implemented  was  linear  discontinuous  (LD). 

This  first  spatial  moment  method  was  chosen  to  show  the  implementation  of 
other  higher  order  spatial  methods  with  the  same  algorithm  as  LC.  Boergers  et 
al  (4:  289-290)  provided  the  angular  flux  representation  for  the  LD  equations.  A 
derivation  for  the  LD  cell  current  equations  for  ordinate  n  with  7Jn  >  /Jn  >  0  and 
ocn<  1 ,  is  presented  in  Appendix  D. 

The  following  definitions  are  used  for  the  LD  quadrature: 


,  3  3  an2 

Clyj  —  l  H"  OCyj  +  £v  +  +  , 

*  4  +  ^„  1  +  30^+^ 

(4.139) 

K  =4  +  eyn  > 

(4.140) 

and 

cn  =  1  +  3 an  +  £yn  ■ 

(4.141) 

Again,  the  process  of  determining  the  equations  for  the  k  vectors  used  to  form 
the  diagonal  matrices  K  =  D(T)  is  analogous  to  the  procedure  used  for  SC.  The 
k  vectors  for  the  top  outward  angular  currents  are: 


(9  +  6bn  -  3anbn  +bn2) 

(4.142) 

«A2 

(3  +  bn)(3  +  cn)an  \rjn\ 

(4.143) 
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(A=-a"(b:+1), 

arfincn 

(4.144) 

( rTL  \  _  (-3  +  (a„-l)h„)  ^„ 

\^O0)n  2 1  l  ’ 

^rrn  \Mn 

(4.145) 

(^)„=(^,)„=a>,<V3)- 

anbn 

(4.146) 

nJ  ^  _  4va„(h„+3) 

\^OSX  )n  ,  ? 

anrncn 

(4.147) 

and 

(1T  ,  _  Av(-3  +  (fli;-l)h„) 

\KOSY )n  ~  ,  i 

anbn 

(4.148) 

The  k 

vectors  for  the  right  outward  angular  currents  are: 

n,RB,  _  (3  +  bn )(cn  +  ^an ) \Mn 

\^OI  )n  /  1  1  ’ 

anbncn\nn\ 

(4.149) 

rhRL^  ocn  (c„-  +9  afI+  cn  (3  -  3  au  +  3  an )) 

\KOI  'n  ~  2  ’ 

(4.150) 

tirRB\  _  \Mn  \ 

\^O0  )n  ~  2 1  1  ’ 

ancn  yln\ 

(4.151) 

(C>»=-(c”t3“"), 

arrncn 

(4.152) 

n,R  ^  _fUR  ^  _  Ay(cn+2an)\Vn\ 

\^OSA)n  \^OEA)n  \  \  ’ 

A 

(4.153) 

,  Av(fl„c„  -  (cn  +  3an ))  //„ 

\KOSX  )n  2 1  l  5 

yln\ 

(4.154) 

and 

nR  ,  Av(c;i+3a„)|//„| 

ykOSY>n  ~  ,ii- 

Kl 

(4.155) 

The  k 

vectors  for  the  top  outward  first  moment  of  the  angular  currents 

are: 
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and 


nJBx  _3{3  +  bn)an 

\^6I  )n  ,  5 

anrncn 

(4.156) 

_  3an  (~ancn  +  (3  +  cn  )an )  \r/n  \ 

>n  "’ll  ’ 

ancn  \Mn\ 

(4.157) 

^ j^TB  ^  _  iP'tfin  3 0Cn  ) 

ancn 

(4.158) 

(*g).~  3“"l'l"l|. 

^iPnPn  \Mn  \ 

(4.159) 

(*L)»  =  (*««)»  =J^L. 

ancn 

(4.160) 

{]t  )  _Ay(ancn-3an2) 

ancn 

(4.161) 

(kT  x  _  Ay3an 

\k6SY)n  ~  , 

a„b„c„ 

(4.162) 

The  k  vectors  for  the  right  outward  first  moment  of  the  angular  currents  are: 


r,RB^  _  3(-3  +  (o„  -l)h„)  nn\ 

yK6I  hr  ,  2 1  l  ’ 

aA  m 

(4.163) 

si  RL\  _3a„(3  +  c„) 

\Kei  )n  7  5 

anPncn 

(4.164) 

ruRB\  _  3or„|//„| 

\Kee  )n  ,ii? 

anbncn  \%\ 

(4.165) 

,uRL\  _(~3  +  «A) 

y^OO  )n  2  5 

anhn 

(4.166) 

(*L),=(0»=^f 4. 

«A  Al 

(4.167) 

,  3Ayan  \fJ.„ 

y^osx'n  ,  I  ’ 

«AC«  *7„| 

(4.168) 
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and 


nR  \  _A y(anbn~3)\Vn 

\K6SY)n  ~  Th 


The  k  vectors  for  the  average  angular  fluxes  are: 

(uB  x  _  (3  +  h„) 

\KAl)n  -  ;  ~  > 

®rfin  ^ln 


(d;),  =  g"<3  +  ‘") 


I /“h  I 


and 


(d«)„=- 


«»)„=- 


(^asa)h  = (^aea) 
(^ASx)n  =~ 
(^ASy)h  =  ~ 


yin  \ 

1 

«Akl  ’ 

Av 


Avar,, 
|^7// 1 

Av 

Vn 


The  k  vectors  for  the  x-moment  of  the  angular  fluxes  are: 

nrB  \  _  3(3  +  bn)cxn 

\KXI  hi  ~  l  “  ) 

®nrrPn  ^ln 


(hL,  _3an(-ancn+(3  +  cn)an) 

\kXl)n  ~  2 

Mn 


(k-xd )  =  (a"C"  ~>^n  ^ 
ancn  Vn 


(k-xo )  =  ^a"c"  ~>(*n  ^ 

ancn  ^ln 


(4.169) 

(4.170) 

(4.171) 

(4.172) 

(4.173) 

(4.174) 

(4.175) 

(4.176) 

(4.177) 

(4.178) 

(4.179) 

(4.180) 
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(k XSaXi  ~  (k x E a) n  ~  '  ||) 

(4.181) 

yin  \ 

n,  ^  _  Ay(a„c„-3a„2) 

y^XSX  )n  2 1  l  ’ 

(4.182) 

ancn  A| 

_  Ay3  a„ 

\KXSY>n  ~  ,  „  • 

(4.183) 

^nbn^n  'Hn 

and 

The  k  vectors  for  the  y-moment  of  the  angular  fluxes  are: 

,_(9  +  3Z>„-3flA) 


and 


kyj  (n,  n)  =  ■ 


k'Y, (n, n)  ■ 


kY0(n,n)  =  - 


a  b~  ri 

n  n  hi 


3a,X3+c„) , 

anbnCn\Vn\  ’ 

3  a„ 


kYe(n,n )  = 


_(-3  +  anbn) 


( bYSA)n  ~  ( bYEA)n  ~ 


i  2  1  I  ^ 

«A  W 

Ay3 


Ai  | 


(bYSx)n  ~ 
(k YSY ) 


Ay3a„ 


raartn  -  ,  i  i  > 

®nbn^n  yin  \ 

A v(gA  -3) 
Vn 


(4.184) 

(4.185) 

(4.186) 

(4.187) 

(4.188) 

(4.189) 

(4.190) 


As  with  WDD,  LD  is  treated  like  LC:  each  ordinate  is  evaluated  to 
determine  the  outgoing  face  and  the  respective  transport  coefficient.  For 
ordinates  that  are  not  exiting  the  top  of  the  cell,  the  same  basic  relations  are 
used  with  an  x-y  exchange,  a  right-left,  reversal  or  a  top- bottom  reversal  where 
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appropriate.  Again,  for  the  DI  method,  these  LD  matrices  (or  any  other  linear, 
first  spatial  method)  use  the  same  solver  algorithm  as  LC. 

C.  Partial  Current  Problem 

1.  Zeroth  Spatial  Moment  Partial  Current  Problem 
The  partial  current  problem  is  set  up  to  establish  the  proper  scale  of 
values  across  the  problem.  Setting  up  the  partial  current  problem  with  the 
improved  cell  shape  information  requires  an  array  similar  to  equation  (3.88).  The 
edge  distribution,  C,  ,  is  defined  as  the  current  along  an  ordinate  divided  by  the 
partial  current  found  by  integrating  over  the  ordinates  in  that  direction  or: 


b  R 


Jr 


(4.191) 


I  >//;■  ’ 

n'eR 

where  wn  is  the  angular  quadrature  weight.  Similarly,  an  edge  distribution  can 


be  defined  for  the  top  edge  as  well: 


br 


Jt 


(4.192) 


2>-y ' 

n'eT 

The  left  and  bottom  edges  are  defined  in  the  same  way. 

The  edge  distributions  allow  the  cell  current  shape  information  to  be 
retained  while  solving  the  partial  current  problem.  Using  the  relations  in 
equation  (4.17),  the  edge  current  along  the  right  edge  in  terms  of  the  incoming 
currents  and  emissions  in  a  cell  is: 


Jout  mOlJln  +  mOI  J In  +  mOI  J  In  +  mOI  Jin  mOEA  A 


A  ' 


(4.193) 
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As  was  done  in  slab  geometry,  the  right  edge  current  relation  can  be  transformed 


JL  =  MoiJi + KUi  +  KJJl  +  +  mroeaea  . 


into  an  equivalent  relation  for  the  partial  currents  and  emissions  in  a  cell  as: 

(4.194) 

In  this  case,  the  coefficients  on  the  partial  currents  and  emissions  are  collapsed 
single  values  determined  using  the  quadrature  weights,  the  coefficient  matrices 
and  the  edge  distributions  as  follows: 


r  RL 


rRT 


and 


KeaEa 


■  X  «4 

ne  out 

mgf f)  . 

n 

(4.195) 

:  X  w» 

neout 

1  rl  y-L  \ 

{mOlC  j  > 

n 

(4.196) 

'■  X  Wh 

neout 

mo/?7)  , 

n 

(4.197) 

:  X  wb  ( 

neout 

m %?)  , 

n 

(4.198) 

=  X  A 

<i{^OEaEa)  ■ 

(4.199) 

neout 


Equations  (4.195)  through  (4.198)  can  be  applied  in  a  similar  manner  to 
determine  the  remaining  partial  currents.  Applying  this  to  the  system  of 


equations  described  in  equation  (4.17)  yields: 


JL 

u  Out 

MLB~ 

J1' 

J  In 

1 

bq 

-4  O 

1 _ 

JR 

u  Out 

M% 

Moi 

MR0B 

JR 

J  In 

+ 

MoeaEa 

JT 

u  Out 

Moi 

M™ 

MT0] 

MTB 

JT 

J  In 

MqeaEa 

I 

MBB 

MB0f 

MBo> 

MBB 

JB 

\_J  In  J 

MboeaEa_ 

(4.200) 
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The  partial  currents  in  equation  (4.200)  are  the  outgoing  partial  currents  for  a 
cell.  The  collapsed  matrices  form  coefficients  for  each  cell  in  the  spatial  mesh 
that  will  be  used  in  the  partial  current  problem. 


CCB 

Figure  4.2.  Setting  up  the  partial  current  problem  for  a  two  cell  by  two 
cell  problem  showing  the  ordering  of  the  partial  currents. 


The  ordering  of  the  partial  currents  in  the  partial  current  problem  is 
important  in  keeping  the  problem  manageable.  The  scattering  contribution  from 
the  orthogonal  directions  increases  the  bandwidth  of  the  sparse  matrix.  To  keep 
the  matrix  bandwidth  manageable  and  provide  a  consistent  pattern  to  implement 
into  the  code,  the  partial  current  problem  was  set  up  using  the  ordering  shown  in 
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figure  4.2.  The  currents  are  numbered  across  the  rows  left  to  right  and  up  then 
up  the  columns  from  bottom  to  top  and  right.  The  quantities,  aL  ,  aR,  aT  ,  and 
aB  represent  the  boundary  conditions  on  the  left,  right,  top  and  bottom  side  of 
the  spatial  mesh  respectively.  The  same  boundary  condition  was  applied  to  all 
the  cells  on  a  respective  edge.  Two  boundary  conditions  were  used:  vacuum 
boundaries  (i.e. aR  =0);  and  reflective  boundaries  (i.e. aR  =1).  The  cell  partial 
currents  in  equation  (4.200)  are  rearranged  across  all  the  cells  to  create  a  matrix 
equation  of  the  form  Ax  =  b  by  using  the  relation  between  the  cell  inward  and 


outward  partial  currents.  For  the  problem  shown  in  figure  4.2,  the  matrix  A  is: 
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where,  for  example,  M[4  is  the  coefficient  that  corresponds  to  the  cell  top 
incoming  partial  current  that  contributes  to  the  outgoing  partial  current  J14  . 
These  values  come  directly  from  the  cell  partial  current  equation  shown  in 
equation  (4.200)  and  are  the  collapsed  matrix  coefficients  for  each  cell.  The 
unknowns  x  in  the  relation  are  the  partial  currents  in  the  ordering  shown  in 

figure  4.3:  x  =(j*  J2  J3  ...  .  J22  J23  J24)  and  the  forcing  vector  b , 

is  derived  from  the  cell  emissions  in  the  same  ordering.  The  elements  of  the 
forcing  vector  are  considered  known  values  for  this  problem. 

As  can  be  seen  by  the  matrix,  this  is  a  sparse  matrix  problem  which  grows 
quickly  as  the  number  of  cells  in  the  problem  increase.  To  solve  the  sparse 
matrix  problem,  a  Compaq  Extended  Math  Library  (CXML)  (6:  11-1)  direct 
sparse  matrix  solver  (cxml_dss.f90)  was  used.  Fortunately,  the  library  routine 
did  not  require  actually  creating  the  matrix  explicitly;  data  was  entered  as 
vectors  which  greatly  increased  to  size  of  the  problems  that  could  be  solved. 

Reapportioning  Partial  Current  from  the  Direct  Solver 

The  partial  current  problem  solution,  J PCP ,  from  the  library  routine  can 
then  be  distributed  back  to  the  cell  edge  currents  using  the  original  £  ,  or  edge 
distribution.  This  forms  the  basis  of  the  iterative  method.  With  the  correct  C, 
value,  the  partial  current  solution  does  not  change  from  the  initial  partial  current 
values.  Since  the  correct  £  value  is  not  known  initially,  an  iteration  with  among 
cell  calculations  on  the  cell  edge  values  must  be  used  to  improve  the  current 
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estimate  of  C,  .  To  distribute  the  cell  partial  current  solution  back  to  the  cell 
edge  values,  the  following  relations  are  used  for  the  zeroth  spatial  moment 
methods: 


1L  -  JL  7l 

Jnew  ^  PCP?  5 

(4.201) 

7 r  -  tr  yR 

Jnew  J  PCP?  ? 

(4.202) 

~.T  _  jT  7 T 

Jnew  PCP  r  > 

(4.203) 

and 

7  B  _  rB  ~?B 

Jnew  PCP  b  • 

(4.204) 

2.  First  Spatial  Moment  Partial  Current  Problem 

First  spatial  moment  methods  must  be  handled  differently  due  to  the  first 
spatial  moment  of  the  edge  current  0  that  is  used  for  these  methods.  The 
solution  can  be  found  through  either  solving  two  simultaneous  systems  of 
equations  or  transforming  the  partial  current  system  of  equations  to  eliminate  the 
0  values.  The  second  choice  was  used  in  order  to  allow  the  use  of  the  same 
routine  for  the  partial  current  problem  that  was  used  with  the  zeroth  spatial 
moment  methods.  To  do  this,  a  new  parameter  is  defined: 

A=-.  (4.205) 

Ji 

where  p  is  a  cell  edge  array  containing  the  number  of  ordinates  in  the  angular 
quadrature  set.  This  new  parameter  permits  the  first  spatial  moment  methods 
current  to  be  written  in  a  form  similar  to  equation  (4.193).  In  this  case,  the 
equation  may  be  written  as: 
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(4.206) 


Ton,  =  (m0I  +  mO0  D(pJ)jL  +  (m0/  +  mol  D(pl,))jln  + 

«t  +<e  vCp\MnH<:+<e  D  (PtM„  +<EJA, 

where  D  is  an  operator  that  creates  a  diagonal  matrix  from  a  vector.  A  similar 
procedure  is  done  for  the  remaining  edge  currents.  This  allows  the  system  of 
equations  in  equation  (4.89)  to  be  written: 
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(4.207) 


The  quantity  in  the  parenthesis  can  be  combined  to  form  a  single  matrix: 


7  L 

JOut 

1R 

JOut 

1T 

JOut 

-■B 

_J Out  _ 

^LL 

—  LR 

—  LT 

—  LB 

m  oi 

mo/ 

m  oi 

m  oi 

J  In 

—  RL 

~RR 

—  RT 

—  RB 

~r  j? 

m  oi 

mo/ 

mo/ 

mo/ 

Jin 

+ 

—  TL 

~TR 

~TT 

—  TB 

1T 

m  oi 

mo/ 

mo/ 

mo/ 

J  In 

—  BL 

—  BR 

~BT 

—  BB 

t 

_ i 

m  oi 

m  oi 

mo/ 

mo/ 

...  L  EA 

mOEAh 

mR  FA 
mOEAh 

mT  FA 
mOEAh 

mB  FA 
mOEAh 


(4.208) 


Here  the  m  indicates  the  quantities  in  the  parentheses  for  equation  (4.207).  Now 
equation  (4.208)  is  in  the  same  form  as  equation  (4.17)  and  the  collapsing  for  the 
partial  current  problem  is  done  the  same  as  for  the  zeroth  spatial  moment 


method.  The  first  moment  partial  current  problem  is  identical  to  the  zeroth 
spatial  moment  problem.  Also,  the  cell  partial  current  solution  is  distributed  to 


the  cell  edge  currents  as  shown  in  equations  (4.201)  through  (4.204).  One 
difference  is  the  cell  edge  current  first  spatial  moment  values  0 .  are  distributed 
as  follows: 
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@new  ~  J  pep  D(/2  )C  > 

(4.209) 

~tR  d  — *R  ~z,R 

&new  ~  J pep  )  C,  , 

(4.210) 

6 new  —  J PCP  D(/?  )  C  i 

(4.211) 

6  new  —  J pep  D(pB)CB.  (4.212) 

Note  the  partial  current  problem  does  not  adjust  the  edge  distributions  t,  .  This 
is  done  during  the  among  cell  calculations  using  the  local  cell  coupling  relations. 
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V.  Validation  and  Performance 


The  code  must  be  validated  before  any  comparisons  to  other  methods  can 
be  made.  The  test  plan  was  implemented  in  three  phases:  initial  checks; 
consistency  checks;  and  accuracy  checks.  These  checks  and  their  results  are 
described  in  this  chapter. 

A.  Validation 

1.  Initial  Checks 

Two  key  areas  for  initial  checks  were  for  both  the  spatial  method  and  the 
angular  quadrature.  For  the  angular  quadrature,  the  weights  and  direction 
cosines  were  tested  using  Mathematica  to  compare  the  ability  of  the  angular 
quadrature  to  exactly  integrate  the  functions  1 ,  /J  ,  [T ,  Ll' .  jll4 ,  r] ,  iy  .  rt f  and 
ij4  over  the  interval  -1  to  1. 

Numerical  testing  confirmed  that  cell  balance  equations  were  satisfied  by 
each  spatial  quadrature  method  as  implemented.  Most  errors  would  show  up  as 
violations  of  the  balance  equations  (13:  176-177).  For  the  zeroth  spatial  moment 
methods,  the  particle  balance  equation  for  a  cell  is: 

Ur  ~  Jl  )4 y  +  UT  ~  Jb  )Ax  +  o-AxAyy//(  =  S4AxAy .  (5.1) 

For  the  first  spatial  moment  methods,  additional  balance  equations  were  used. 
The  x  moment  balance  equation: 

3Ur  +  JL  ~ 2M¥a)a. v  +  (0T  ~ &b)Ax  +  o-AxAyy/x  =  SxAxAy  ,  (5.2) 
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and  the  y  moment  balance  equation: 


3 (jT  +  jB  -  2r/y/A)Ax  +  (0R  -  9 L  )Av  +  <7AxAyy/Y  =  SYA\Ay  .  (5-3) 

Both  the  angular  quadrature  testing  and  all  cell  balance  relations  for  all  spatial 
methods  were  confirmed. 

2.  Consistency  Checks 

The  consistency  testing  was  broken  into  two  portions:  symmetry  tests  and 
aspect  ratio  tests. 

Symmetry  Tests 

In  this  phase,  the  testing  validated  that  boundary  conditions  and  indexing 
were  consistently  implemented.  (This  test  identifies  copy-paste-edit  errors.)  The 
quantities,  aL  ,  aR  .  aT  .  and  aB  are  used  to  specify  the  boundary  conditions  on 
the  left,  right,  top  and  bottom  side  of  the  spatial  mesh  respectively.  The  same 
boundary  condition  was  applied  to  all  the  cells  on  a  respective  edge.  Again,  two 
boundary  conditions  were  used:  vacuum  boundaries  (i.e.  aR  =0);  and  reflective 
boundaries  (i.e.  aR  =1).  The  scattering  ratio  is  varied  in  these  tests.  It  is 
defined  as  the  ratio  of  the  scattering  cross  section  ( as )  to  the  total  cross  section 

(Oi): 


as 

c  =  — — 


(5.4) 


Various  symmetries  are  compared  to  ensure  the  same  result  is  calculated  when 


only  the  orientation  of  the  problem  is  changed.  For  the  different  problems 
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examined,  the  cell  average  scalar  flux  should  be  the  same  value  and  the  rate  of 
convergence  should  be  identical. 

The  first  test  problem  in  this  phase  is  the  uniform  universe  test. 

Reflective  boundaries  are  set  on  all  boundaries  and  the  values  are  set  as  shown  in 
figure  5.1. 


aT  =  1 


Figure  5.1.  Problem  values  for  the  uniform  universe  test  problem. 

This  test  problem  was  chosen  because  it  is  one  of  the  few  transport 
problems  with  a  closed  form  solution.  A  flux  solution  is  found  by  integrating  the 
BTE  over  all  angles: 

f  [RV^+(7^  =  (7!(l)  +  £,]fifQ,  (5.5) 

\/a 

but  Vl//  =  0  and  \jf  is  independent  of  Q  for  this  uniform  problem.  This  also 
means  that: 

=  (5-6) 

where  a  —  1  based  on  the  normalization  for  the  angular  quadrature  set  where: 
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(5.7) 


f  dQ  =  1. 
va 


Equation  (5.5)  yields: 


a,(P  =  axli+E . 


(5.8) 


or, 


in,  -  (J.}0  =  fj.,0  =  e 


(5.9) 


where  <Ja  is  the  absorption  cross  section.  The  solution  for  the  scalar  flux  is: 

</>  =  —  .  (5.10) 

Also,  for  the  angular  quadrature  set  normalization,  the  value  of  the 
converged  angular  and  scalar  flux  in  a  cell  should  be  the  same  and  equal  to  — . 

The  next  symmetry  test  problem  examined  the  effect  of  boundary 
conditions  on  the  solution  by  setting  three  sides  of  the  problem  with  reflective 
boundaries  and  the  remaining  side  with  a  vacuum  boundary.  The  side  with  the 
vacuum  boundary  is  rotated  through  all  possible  cases,  and  the  (rotated  or 
reflected)  converged  solutions  should  be  identical  in  each  case.  The  problem 
values  for  this  test  are  shown  in  figure  5.2. 
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aT  =  1 


Figure  5.2.  Problem  values  for  the  single  vacuum  boundary  and  three 

reflective  boundary  problem. 

The  last  symmetry  test  problem  in  this  phase  examines  additional 
rotational  symmetries.  In  this  problem,  two  adjacent  boundaries  are  reflective 
and  the  other  two  are  vacuum,  then  the  boundary  conditions  are  reversed. 

Again,  for  both  of  these  cases,  the  converged  results  should  be  identical.  The 
problem  values  were  the  same  as  figure  5.2  with  the  exception  of  the  vacuum 
boundaries. 

The  results  of  some  of  the  validation  tests  follow.  For  brevity,  the  results 
of  the  WDD  method  are  shown  for  some  of  the  tests,  the  other  spatial 
quadratures  had  similar  results. 

Symmetry  Test  Results 

The  results  of  the  uniform  universe  test  are  shown  in  table  5.1.  The  test 
was  done  using  the  S6  angular  quadrature  and  a  tolerance  of  10-5 .  Each  spatial 
method  converged  in  one  iteration  for  the  DI  method.  As  noted  earlier,  the  value 
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of  the  scalar  flux  in  a  cell  can  be  calculated  and  for  this  test  the  value  should  be 


0.9325873.  An  independent  source  iteration  (SI)  solution  was  also  done  for 
comparison.  Note  that  while  the  SI  solution  meets  the  requested  tolerance,  it 
does  not  have  the  precision  that  the  DI  methods  have  for  this  solution.  The 
number  of  SI  iterations  required  to  meet  the  same  tolerance  is  listed,  which  is 
significant  for  a  relatively  simple  test  problem.  In  addition,  for  this  test,  the 
average  angular  flux  in  a  cell  had  the  same  value  as  the  scalar  flux  as  expected. 


Table  5.1.  Results  of  the  uniform  universe  test. 


Spatial 

Distribution  Iteration 

Source  Iteration 

method 

Scalar  Flux 

Number  of 

Scalar  Flux 

Number  of 

iterations 

iterations 

WDD 

0.9325873 

1 

0.9325859 

38 

SC 

0.9325873 

1 

0.9325859 

35 

LD 

0.9325873 

1 

0.9325871 

21 

LC 

0.9325873 

1 

0.9325871 

21 

The  results  of  the  one  vacuum  boundary  symmetry  tests  showed  that  the 
problems  returned  identical  values  for  the  scalar  flux,  iterations  to  convergence 
and  maximum  and  minimum  scalar  flux  values  for  each  vacuum  boundary 
location.  The  test  was  also  done  using  the  S6  angular  quadrature  and  a 
tolerance  of  10-5  .  As  noted  earlier,  the  results  of  each  different  vacuum 
boundary  should  be  identical  as  the  vacuum  boundary  is  rotated  around  the 
problem  grid  if  the  boundary  conditions  and  indexing  are  correct  for  either  a 
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right /left  or  t-op/bottom  exchange.  The  other  spatial  quadratures  have  similarly 
identical  results. 

In  addition,  the  average  angular  flux  in  each  cell  was  compared  for  the 
results  for  the  right/left  and  t-op/bottom  tests  respectively  by  exchanging  the 
array  indices  for  the  right/top  and  comparing  this  to  the  left /bottom  test  results. 
The  exchanged  right  vacuum  boundary  cell  average  angular  flux  values  compared 
to  the  test  left  vacuum  boundary  cell  average  angular  flux  with  a  SRD  of 
1.28x10  L’ .  Also,  the  exchanged  top  vacuum  boundary  cell  average  angular  flux 
values  compared  to  the  test  bottom  vacuum  boundary  cell  average  angular  flux 
with  a  SRD  of  1.28x10  13 .  The  other  spatial  quadratures  have  similar  results. 
These  test  results  all  used  a  convergence  tolerance  of  10-5 ,  and  the  SRD  is 
consistent  with  the  rounding  errors  associated  with  the  machine  arithmetic  for 
the  different  test  solutions. 

The  next  symmetry  tests,  two  vacuum  boundaries,  also  returned  identical 
values  for  the  scalar  flux,  iterations  to  convergence  and  maximum  and  minimum 
scalar  flux  values  for  both  vacuum  boundary  cases.  Computations  for  this  test 
used  the  S6  angular  quadrature  and  a  tolerance  of  10-5.  For  this  test,  the 
results  for  the  different  vacuum  boundaries  should  also  have  been  identical  as  the 
vacuum  boundaries  are  rotated  on  the  problem  grid  if  boundary  conditions  and 
indexing  for  an  x/y  exchange  are  properly  implemented.  The  other  spatial 
quadratures  again  had  similar  identical  results. 
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For  this  case  as  well,  the  average  angular  flux  in  each  cell  was  compared 
by  doing  another  exchange  of  the  array  indices  for  the  right/top  test  results. 

The  exchanged  right/top  vacuum  boundary  cell  average  angular  flux  values 
compared  to  the  test  left/bottom  vacuum  boundary  cell  average  angular  flux 
with  a  SRD  of  1.60x10  13 ,  which  is  consistent  with  rounding  errors  for  the 
different  tests.  The  other  spatial  quadratures  had  similar  results. 

Aspect  Ratio  Tests 

The  symmetry  test  problems  used  a  40x40  grid  of  square  cells.  The  next 
series  of  test  problems  in  this  phase  of  testing  uses  various  aspect  ratios  Ay :  Ax 
while  keeping  the  cross  section  and  boundary  conditions  the  same  as  shown  in 
figure  5.3.  Again,  the  results  for  the  converged  solution  should  be  identical  when 
x  and  y  values  are  interchanged.  Aspect  ratios  of  1:2,  1:4  and  1:8  were  compared 
to  aspect  ratios  of  2:1,  4:1  and  8:1. 

The  aspect  ratio  tests  returned  identical  values  for  the  scalar  flux, 
iterations  to  convergence  and  maximum  and  minimum  scalar  flux  values  for  all 
spatial  quadratures.  Again,  the  test  was  also  done  using  the  S6  angular 
quadrature  and  a  tolerance  of  1 0  5 .  As  noted  earlier,  this  test  confirms  that  cells 
with  aspect  ratios  other  than  one  returned  consistent  results  when  x  and  y  values 
are  interchanged. 


102 


3.  Source  Iteration  Comparison 


The  last  series  of  tests  in  this  phase  compare  the  converged  solution  from 
conventional  source  iteration  with  the  DI  solution.  The  scattering  ratio  was  kept 
low  so  that  the  SI  solution  would  not  suffer  from  false  convergence.  An  example 
of  the  problem  and  boundary  conditions  used  is  shown  in  figure  5.3  for  a  10x10 
spatial  mesh.  The  spatial  mesh  is  refined,  from  10  cells  by  10  cells  to  100  cells  by 
100  cells.  The  cell  scalar  flux  results  for  both  SI  and  DI  are  compared  to  ensure 
the  converged  results  are  consistent. 


aT  =  0 


Figure  5.3.  Problem  variables  for  the  source  iteration  comparison  test 
problem. 

Source  Iteration  Test  Results 

The  results  of  the  source  iteration  comparison  tests  are  shown  in  figure 
5.4.  Again,  the  test  was  done  using  the  S6  angular  quadrature  and  a  tolerance  of 
10  5 .  The  solid  line  is  the  requested  tolerance  of  10~5.  As  noted  earlier,  this  test 
confirmed  that  an  independent  source  iteration  calculation  returned  the  same 
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values  for  the  cell  scalar  fluxes  as  the  DI  method.  The  SRD((/>di  ,(/>SI)  was  less 
than  3.6 xlCT6,  which  is  less  than  the  1 0  tolerance,  in  every  cell  for  seven  trials 
with  grids  ranging  from  10x10  to  100x100.  This  gives  confidence  that  the  code  is 
consistent,  however  it  is  still  possible  that  both  the  DI  and  the  SI  codes  could  be 
off  by  a  common  factor.  To  eliminate  this  possibility,  the  results  are  next 
compared  to  an  independent  solution. 

4.  Benchmarking 

After  the  initial  checks  and  consistency  checks  are  done,  it  is  evident  that 
the  results  from  the  code  are  consistent.  Getting  the  same  results  for  different 
problems  from  two  different  methods  within  the  code  shows  consistency,  but  it 
does  not  demonstrate  accuracy.  To  do  this,  the  converged  results  must  be 
compared  to  a  known  solution  (benchmarked).  Mathews’  vacuum  duct  problem 
(10:  x-8)  is  used  as  a  benchmark.  The  benchmark  problem  is  shown  in  figure  5.5. 


<XB=  0 

Figure  5.4.  Problem  variables  for  the  benchmark  problem. 
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Benchmark  Test  Results 


The  results  of  a  benchmark  test  for  the  SC  spatial  quadrature  is  shown  in 
figure  5.5.  Again,  the  test  was  also  done  using  the  S8  angular  quadrature  and  a 
tolerance  of  10-' .  The  solid  line  shows  an  independent  Monte  Carlo  solution  to 
the  same  problem  (10:  x-8).  A  ray  effect  due  to  the  use  of  the  S8  angular 
quadrature  is  seen  in  the  location  of  the  peak  of  the  graph.  Physically,  the  peak 
should  be  located  over  the  duct,  as  shown  by  the  Monte  Carlo  solution.  This  ray 
effect  behavior  is  also  consistent  with  previous  results  for  this  problem  (10:  x- 
10). 
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Figure  5.5.  Results  of  the  SC  comparison  with  a  Monte  Carlo  solution  to 
the  benchmark  problem.  The  plot  of  the  partial  current  out  the  top  edge 
is  shown  for  both  methods. 
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The  results  of  a  benchmark  test  for  linear  discontinuous  is  shown  in  figure 


5.6.  Again,  the  test  was  done  using  the  S$  angular  quadrature  and  a  tolerance  of 
1(T5 .  The  solid  line  shows  an  independent  Monte  Carlo  solution  to  the  same 
problem.  As  with  the  SC  solution,  a  ray  effect  due  to  the  use  of  the  S8  angular 
quadrature  is  seen  in  the  location  of  the  peak  of  the  graph.  Again,  this  ray  effect 
behavior  is  also  consistent  with  previous  results  for  this  problem. 
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Figure  5.6.  Results  of  the  LD  comparison  with  a  Monte  Carlo  solution  to 
the  benchmark  problem.  The  plot  of  the  partial  current  out  the  top  edge 
is  shown  for  both  methods. 

For  both  the  SC  and  LD  solutions  shown  in  figures  5.5  and  5.6,  the 
important  observation  is  not  the  ray  effect,  but  the  magnitude  of  the  partial 
current  calculated  for  both  spatial  methods.  The  scale  of  the  DI  result  is 
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comparable  to  the  Monte  Carlo  result  in  either  case  and  validates  the  accuracy  of 


the  comparison  done  with  SI  in  figure  5.4. 

B.  Routine  Problem  Comparison 

This  section  demonstrates  the  efficiency  of  this  method  by  comparing  the 
results  for  the  DI  method  to  published  results  for  DSA  methods  for  three 
different  problems:  varying  aspect  ratios,  varying  scattering  ratios  and  varying 
mesh  size,  for  given  problems.  The  three  problems  were  not  particularly 
challenging  for  either  method,  but  show  how  DSA  methods  and  SI  methods 
compare  to  the  DI  method  for  relatively  straightforward  problems.  Morel  et  al. 
(14:  309-10)  published  results  for  DSA  using  Bi-Linear  Nodal  (BLN)  and  Waring 
et  al.  (17:  124-25)  published  results  for  DSA  using  Linear  Bi-Linear  Nodal 
(LBLN)  spatial  methods  for  the  same  set  of  three  problems.  The  comparative 
measure  used  for  each  problem  is  the  number  of  iterations  needed  to  converge  the 
cell  scalar  flux  to  a  given  tolerance.  The  DSA  methods  have  an  inner  loop  which 
is  used  to  estimate  the  residual  error  at  each  step.  For  these  problems  listed,  the 
DSA  inner  loop  used  a  minimum  of  three  passes  to  update  the  residual  error 
estimate,  while  the  DI  method  only  had  one  pass  through  the  among  cell 
calculations  (14:  306).  However,  for  comparison  purposes,  an  iteration  is  one 
complete  cycle  in  each  case,  which  should  be  a  conservative  comparison  for  the 
DI  method. 
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1.  Aspect  Ratios 


The  first  comparison  problem  examines  results  for  a  set  of  grids  which 
differ  in  cell  aspect  ratio,  Ay :  Ax  .  The  basic  parameters  for  a  homogeneous 
medium  problem  (14:  309,  17:  124)  are  shown  in  figure  5.7.  The  problem  is 
intended  to  show  for  DSA  methods  the  effectiveness  in  terms  of  error  reduction 
per  iteration.  The  problem  was  done  using  an  S4  angular  quadrature  and 
converged  to  a  tolerance  of  1 0  4  using  the  cell  average  scalar  flux.  The  spatial 
grid  has  25x25  rectangular  cells  in  each  case;  the  problem  size  differs  among  the 
cases.  The  cells  are  not  necessarily  square.  Aspect  ratios  of  Av :  Ax  =  2  : 1 ,  5:1, 
10:1  and  20:1  were  tested. 


aT=  0 


Figure  5.7.  Problem  variables  for  the  DSA  aspect  ratio  comparison. 

The  results  for  the  aspect  ratio  tests  for  both  DI  and  DSA  (14:  309,  17: 
124)  are  shown  in  table  5.2.  As  can  be  seen  in  the  table,  the  zeroth  spatial 
moment  methods  using  DI  (WDD  and  SC)  converged  faster  than  the  DSA 
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methods,  while  the  first  moment  methods  using  DI  (LD  and  LC)  were 
comparable  to  the  DSA  methods.  The  stability  of  the  DSA  methods  shows  that 
the  iteration  count  does  not  increase  as  the  aspect  ratio  increases.  The  zeroth 
spatial  moment  methods  also  show  this,  while  the  DI  method  shows  a  slight 
increase  for  high  aspect  ratios  for  first  moment  methods.  For  the  20:1  case,  DI 
takes  more  iterations  (10  for  LD,  7  for  LC)  than  DSA  (6  for  BLN,  6  for  LBLN) 
but  uses  fewer  discrete  ordinates  sweeps  (10  for  LD,  7  for  LC)  than  DSA(  >  18 
for  each  DSA  calculation).  While  this  problem  does  not  definitely  show  the  DI 
method  as  better,  it  does  show  that  DI  requires  of  the  same  order  of  iterations  to 
converge  for  a  totally  scattering  problem. 


Table  5.2.  DSA  aspect  ratio  comparison  results 


DI  Methods 

DSA 

Methods 

Ax 

LBLN 

1.0 

1.0 

4 

5 

5 

5 

8 

8 

1.0 

5.0 

3 

3 

5 

5 

8 

8 

1.0 

10.0 

3 

3 

8 

8 

8 

8 

5.0 

5.0 

2 

2 

3 

5 

6 

6 

5.0 

10.0 

2 

2 

3 

4 

6 

6 

5.0 

100.0 

2 

2 

10 

7 

6 

6 

10.0 

10.0 

2 

2 

3 

5 

5 

5 

10.0 

100.0 

2 

2 

5 

5 

5 

5 

100.0 

o 

o 

o 

2 

2 

4 

5 

5 

5 
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The  same  test  was  done  using  an  adaptive  DO  sweeping  technique  in 


shown  chapter  two  and  described  later  in  this  chapter.  The  results  are  shown  in 
table  5.3  with  the  original  DSA  results  for  comparison.  The  additional  DO 
sweeps  for  each  iteration  only  decrease  the  iteration  count  for  a  few  of  the  tests, 
but  this  is  by  design.  The  adaptive  technique  is  only  to  use  additional  DO 
sweeps  for  an  iteration  where  the  DI  method  is  converging  relatively  slowly, 
which  is  only  two  of  the  tests  the  1:10  and  1:  20  cases.  For  these  problems,  the 
number  of  iterations  to  convergence  is  a  third  smaller.  For  the  other  cases,  the 
iterations  to  convergence  is  about  the  same  or  one  less. 


Table  5.3.  DSA  aspect  ratio  comparison  adaptive 
DO  sweep  results 


DI  Methods 

DSA 

Methods 

Ax 

LBLN 

1.0 

1.0 

4 

5 

4 

3 

8 

8 

1.0 

5.0 

3 

3 

4 

3 

8 

8 

1.0 

10.0 

3 

3 

5 

5 

8 

8 

5.0 

5.0 

2 

2 

3 

3 

6 

6 

5.0 

10.0 

2 

2 

3 

3 

6 

6 

5.0 

100.0 

2 

2 

4 

3 

6 

6 

10.0 

10.0 

2 

2 

3 

3 

5 

5 

10.0 

100.0 

2 

2 

4 

4 

5 

5 

100.0 

i—1 

o 

o 

o 

2 

2 

2 

3 

5 

5 
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The  adaptive  DO  sweeping  technique  shows  the  DI  method  to  be  slightly 
better  than  DSA  for  this  case  in  terms  of  iterations  to  reach  convergence  for  this 
problem. 

Scattering  Ratios 

The  next  comparison  problem  examines  results  for  grids  which  differ  in 
cell  scattering  ratio.  The  basic  parameters  for  another  homogeneous  medium  are 
shown  (14:  309,  17:  124)  in  figure  5.8.  The  problem  is  intended  to  show,  for  DSA 
and  source  iteration  methods,  the  dependence  of  the  efficiency  upon  the 
scattering  ratio.  The  problems  were  solved  using  an  S4  angular  quadrature  and 
cell  average  scalar  fluxes  were  converged  to  a  tolerance  of  10  4 .  The  problem 
uses  a  25x25  cell  grid  with  Ax  =  Ay  =  1  mean  free  path  (mfp). 


Figure  5.8.  Problem  variables  for  the  DSA  scattering  ratio  comparison. 

The  results  for  the  scattering  ratio  comparison  for  unaccelerated  SI,  for 
DI,  and  for  DSA-SI  are  shown  in  table  5.3.  First  consider  SI  versus  DI.  SC  and 
LC  have  similar  relative  performance  for  SI  and  DI  as  WDD  and  LD.  Both  the 
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zeroth  spatial  moment  methods  using  DI  converged  faster  than  source  iteration 


using  the  same  spatial  method  and  angular  quadrature.  The  DI  zeroth  spatial 
moment  methods  show  that  the  iteration  count  went  from  three  to  four  as  the 
scattering  ratio  increased  for  this  problem,  while  the  SI  methods  climbed  from 
twenty  to  over  two  thousand  with  higher  scattering  ratios.  The  DI  first  spatial 
moment  methods  again  show  the  iteration  count  went  from  three  to  five  as  the 
scattering  ratio  increased  for  this  problem,  while  the  SI  methods  increased  from 
ten  to  over  two  thousand  with  higher  scattering  ratios.  It  also  shows  the 
advantage  of  first  moment  methods  for  source  iteration:  the  iteration  count  is 
much  lower  for  the  same  problem  than  with  a  zeroth  moment  source  iteration 
method.  For  DI  methods,  the  iteration  count  was  almost  identical  for  both 
zeroth  and  first  moment  methods.  This  problem  demonstrates  the  DI  method  as 
superior  to  (unaccelerated)  source  iteration  for  this  case. 
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Table  5.4.  DSA  scattering  ratio  comparison  results. 


SI 

DI 

DSA 

c 

WDD 

LD 

WDD 

LD 

BLN 

LBLN 

1.0 

2379 

2020 

4 

5 

8 

8 

0.9 

171 

94 

3 

4 

7 

7 

0.8 

93 

50 

3 

4 

7 

6 

0.7 

64 

34 

3 

4 

6 

6 

0.6 

49 

26 

3 

4 

5 

5 

0.5 

40 

20 

3 

3 

5 

5 

0.4 

34 

17 

3 

3 

5 

4 

0.3 

29 

14 

3 

3 

4 

4 

0.2 

26 

12 

3 

3 

4 

4 

0.1 

23 

9 

3 

3 

3 

3 

Next  consider  DI  versus  DSA.  The  DI  method  converged  slightly  faster 
than  DSA  for  almost  all  scattering  ratios.  The  DI  methods  show  very  little 
increase  in  iteration  (from  3  to  5)  with  increasing  scattering  ratio  for  this 
problem,  while  the  DSA  methods  show  a  larger  increase  (from  3  to  8).  This 
problem  also  shows  the  DI  method  to  be  slightly  better  than  DSA  for  this  case  in 
terms  of  iterations  to  reach  convergence. 

3.  Two  Material  Problem 

The  last  comparison  problem  examines  results  for  grids  which  differ  in  spatial 
mesh  refinement  for  a  two  material  problem.  The  basic  parameters  for  another 
homogeneous  medium  are  shown  (14:  309,  17:  124)  in  figure  5.9.  The  problem  is 
intended  to  show  the  effectiveness  of  DSA  for  inhomogeneous  problems.  The 
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problems  were  solved  using  an  S4  angular  quadrature  and  cell  average  scalar 
fluxes  were  converged  to  a  tolerance  of  1 0  4 .  The  spatial  grid  size  remained  fixed 
at  50  cm  for  this  problem  while  Ax  and  Ay  both  vary  at  the  same  ratio,  thereby 
refining  the  spatial  mesh  for  the  problem.  Mesh  sizes  of  5x5,  10x10,  25x25  and 
50x50  were  tested.  For  these  mesh  sizes,  the  cell  thicknesses  were  10  cm,  5  cm,  2 
cm  and  1  cm  respectively. 


aT=  0 


ccB  —  1 

Figure  5.9.  Problem  variables  for  the  two  material  DSA  comparison 
problem. 

The  DI  results  for  the  mesh  refinement  problem,  along  with  the  published 
DSA  results  (14:  309,  17:  124),  are  shown  in  table  5.5.  As  can  be  seen  in  the 
table,  the  DI  method  converged  slightly  faster  than  the  DSA  methods  for  these 
cell  sizes.  Again,  while  this  problem  does  not  definitely  show  the  DI  method  as 
better  than  DSA,  it  does  again  show  that  DI  converges  in  the  same  number  of 
iteration  or  slightly  fewer  iterations  for  a  highly  scattering  problem. 
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Table  5.5.  DSA  two  material  comparison  results. 


DI  Methods 

DSA  Methods 

Mesh 

WDD 

SC 

LD 

LC 

BLD 

LBLN 

5x5 

2 

2  I 

4 

5 

6 

6 

10x10 

2 

3 

4  1 

5 

8 

7 

25x25 

3 

4  1 

4  ! 

4 

9 

8 

50x50 

4 

6  1 

6  | 

5 

8 

7 

In  chapter  two,  figure  2.1  shows  an  inner  loop  doing  the  local  balance 
coupling  labeled  “iterate  as  needed”.  The  discussion  following  the  figure 
discussed  the  fact  that  some  problems  needed  additional  loops  with  discrete 
ordinates  sweeping.  For  the  problems  presented  so  far,  only  one  discrete 
ordinates  sweep  was  sufficient  for  the  DI  method  to  converge  in  a  few  iterations. 
However,  there  were  problems  in  which  additional  loops  through  the  discrete 
ordinates  sweeping  were  needed  but  this  also  depended  on  the  spatial  method 
used.  This  led  to  an  adaptive  technique  which  varied  between  one  and  ten 
sweeps  depending  on  the  properties  of  the  problem.  Timing  analysis  showed  that 
ten  sweeps  would  at  most  double  the  time  for  an  iteration.  For  each  sweep,  the 
scattering  source  was  updated  using  the  cell  edge  currents  and  the  scattering 
source  was  used  to  calculate  new  cell  edge  currents.  A  detailed  analysis  of  the 
adaptive  technique  is  presented  in  chapter  seven,  but  the  technique  was  used  for 
some  of  the  problems  in  chapter  six,  as  well  as  in  the  first  DSA  comparison 
problem  for  this  chapter. 
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This  chapter  showed  the  validity  of  the  DI  results  through  a  variety  of 


test  problems.  In  addition,  this  chapter  showed  that  the  DI  method  performed 
much  better  than  SI  for  higher  of  scattering  ratios.  Finally,  these  problems  show 
the  performance  of  DI  is  comparable  to  the  effectiveness  of  DSA  with  a  similar 
computational  effort  based  on  a  conservative  iteration  count.  In  the  next- 
chapter,  problems  where  synthetic  acceleration  has  difficulties  are  examined. 
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VI.  Challenging  Problems  -  Comparison  with  DSA  and  TSA 


A.  Where  DSA  Loses  Effectiveness 

The  previous  chapter  shows  how  DI  performance  was  comparable  to  DSA 
on  several  routine  problems.  Recently,  it  has  been  have  shown  that-  DSA  can 
lose  its  effectiveness  or  converges  slowly  for  a  particular  problem. (3:  213,  18:  1) 
This  was  shown  using  a  test  that-  has  alternating  layers  of  two  different  materials 
that  are  highly  scattering.  The  particular  problem’s  parameters  given  by  Azmy 
(3:  228-229)  are  shown  in  figure  6.1.  For  this  problem,  different  total  cross 
sections  will  be  compared  with  different  mesh  sizes  which  varied  from  a  10x10 
spatial  mesh  to  a  160x160  spatial  mesh. 

aT  =  0 


c  =  l 


aR  =  0 

□  - 

□  y* 

< —  Ax — ►  ^  y<* 

Source  Region 

Figure  6.1.  Problem  variables  for  the  Azmy  Periodic  Horizontal  Interface 
(PHI). 
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An  S6  angular  quadrature  was  used  for  this  problem  to  compare  with  the 
published  results  in  the  article.  The  convergence  tolerance  of  10  6  was  used  for 
the  relative  difference  in  the  cell  average  scalar  flux. 

The  measure  of  effectiveness  used  for  this  problem  is  the  spectral  radius. 
For  a  converging  system  of  equations,  the  spectral  radius  is  between  zero  and 
one.  A  spectral  radius  which  is  close  to  zero  indicates  that  the  system  of 
equations  converges  rapidly.  Conversely,  a  spectral  radius  close  to  one  indicates 
the  system  of  equations  converges  slowly.  Additionally,  a  spectral  radius  of  one 
or  greater  indicates  system  of  equations  that  diverges  (6:  229).  Often  calculating 
the  eigenvalues  or  spectral  radius  for  a  large  system  of  equations  is  impractical. 
Azmy  estimates  the  spectral  radius  using  the  ratio  of  the  L2  norm  of  the  residual 
in  the  cell  average  scalar  flux  to  the  previous  iterate  as  follows: 

0/-1  —  2 1|2 

The  spectral  radius  is  computed  for  the  iteration  in  which  the  problem  met  the 
convergence  tolerance  (3:  213-216). 

For  the  DI  method,  testing  showed  that  the  spectral  radius  calculated  this 
way  could  vary  with  the  chosen  tolerance  or  iteration  even  though  the  method 
was  converging  in  a  few  iterations.  Another  method  of  estimating  the  spectral 
radius  or  convergence  rate  was  developed.  The  maximum  SRD  of  the  scalar  flux 


(6.1) 
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between  iterations  is  shown  for  two  different  cross  sections  in  figures  6.2  and  6.3. 


The  solid  line  represents  the  convergence  tolerance  in  Azmy’s  problem. 
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Figure  6.2.  Convergence  rates  for  the  Azmy  Periodic  Horizontal  Interface 
(PHI)  problem  with  DI  using  WDD  and  S6  at  <7  =  10  cross  section  for 
various  mesh  sizes. 


These  problems  are  done  to  a  much  tighter  tolerance,  10.  14  and  for  all  the 
different  spatial  meshes  that  the  DSA  test  problem  were  done.  The  figures 
demonstrate  two  points,  the  maximum  SRD  of  the  scalar  flux  decreases  by  a 
fairly  constant  amount  per  iteration  and  the  problem  can  be  run  to  very  tight 
tolerances  which  show  the  problem  does  not  suffer  from  bad  numerical 
conditioning.  The  rate  of  decrease  in  the  maximum  SRD  is  the  DI  method 
estimate  of  the  spectral  radius  or  convergence  rate.  The  convergence  rate  is 
found  by  doing  a  linear  regression  of  the  linearized  data  which  is  shown  in  the 
figures.  One  note  is  that  this  maximum  SRD  estimate  is  an  asymptotic  value. 


119 


For  problems  with  reasonable  tolerances,  for  example  10  4  is  commonly  used,  the 


problem  would  converge  faster  than  the  DI  estimate  of  the  spectral  radius  would 


predict. 
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Figure  6.3.  Convergence  rates  for  the  Azmy  Periodic  Horizontal  Interface 
(PHI)  problem  with  DI  using  WDD  and  S6  cross  section  &  -  160  for 
various  mesh  sizes. 


Another  observation  from  the  two  plots  shown  in  figure  6.2  and  6.3  is  that 
the  rate  of  convergence  does  change  with  cross  sections  and  mesh  size.  For  the 
cross  section  shown  in  figure  6.2  the  spectral  radius  is  fairly  constant  as  the  mesh 
is  refined.  This  can  be  seen  by  the  fact  the  iteration  count  does  not  change 
significantly  as  the  spatial  mesh  is  changed.  On  the  other  hand,  for  the  cross 
section  used  in  figure  6.3,  the  number  of  iterations  needed  to  reach  the  final 
tolerance  almost  doubles  as  the  mesh  gets  larger.  The  convergence  rates  for  these 
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plots  will  be  discussed  shortly.  For  the  DI  results  shown  in  table  6.2  and  6.4,  the 


adaptive  DO  sweeping  technique  was  used. 

1.  Weighted  Diamond  Difference  Comparison 

The  published  results  for  the  Azmy  PHI  problem  using  DSA  with  WDD 
and  an  S6  angular  quadrature  are  shown  in  table  6.1  (3:  231).  As  seen  in  the 
table,  the  spectral  radius  increases  strongly  with  the  number  of  cells,  indicating 
slower  convergence.  For  this  DSA  method,  going  to  larger  problems  of  this  type 
will  lead  to  slowly  converging  solutions.  Hence,  DSA  is  no  longer  accelerating 
the  solution  for  a  large  enough  problem  of  this  type  or  loses  effectiveness. 

Table  6.1.  Published  DSA  with  WDD  Results. 


Cross  Sections 

Mesh 

10 

20 

40 

80 

160 

10x10 

0.100 

0.039 

0.010 

0.002 

4E-4 

20x20 

0.241 

0.132 

0.044 

0.010 

0.002 

40x40 

0.422 

0.316 

0.151 

0.046 

0.010 

80x80 

0.581 

0.539 

0.360 

0.160 

0.048 

160x160 

0.683 

0.713 

0.609 

0.386 

0.165 

The  DI  results  for  the  Azmy  PHI  problem  using  WDD  with  DI  and  S6 
angular  quadrature  are  shown  in  table  6.2.  The  spectral  radii,  or  convergence 
rates,  listed  in  the  table  were  determined  using  the  slope  of  the  linearized  plots, 
as  described  previously.  Contrary  to  the  DSA  solutions,  the  DI  method  does  not 
increase  strongly  for  larger  problems.  In  addition,  the  total  number  of  iterations 
needed  to  solve  a  difficult  problem  remains  small.  For  this  problem,  DI  with 
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WDD  showed  good  performance  and  was  considerably  better  than  DSA  with  the 


same  spatial  and  angular  quadratures. 

Table  6.2.  DI  with  WDD  results. 


Cross  Section 

Mesh 

10 

20 

40 

80 

160 

10x10 

0.079671 

20x20 

0.077822 

0.016199 

0.002528 

40x40 

0.09177 

80x80 

0.08531 

0.081133 

0.036083 

160x160 

0.122462 

0.064091 

0.081433 

2.  Linear  Discontinuous  Comparison 

Azmy’s  results  for  the  Azmy  PHI  problem  using  a  Bi-Linear  Nodal  method 
with  DSA  and  S6  angular  quadrature  are  shown  in  table  6.3  (3:  232).  Again,  the 
spectral  radius  increases  for  larger  meshes  for  certain  cross  sections.  For  this 
DSA  method,  going  to  larger  problems  will  lead  to  slowly  converging  solutions. 
For  example,  the  spectral  radius  listed  for  the  cross  section  (7  =  20  and  mesh  of 
160,  took  388  iterations  to  meet  the  tolerance  of  10~6  (3:  231).  In  addition,  for 
several  cross  sections,  this  method  diverged. 
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Table  6.3.  Published  DSA  with  BLN  Results. 


Cross  Sections 

Mesh 

10 

20 

40 

80 

160 

10x10 

0.355 

0.254 

0.192 

D 

D 

20x20 

0.543 

0.417 

0.317 

D 

D 

40x40 

0.717 

0.607 

0.452 

D 

D 

80x80 

0.836 

0.688 

0.624 

D 

D 

160x160 

0.901 

0.392 

0.671 

D 

D 

The  DI  results  for  the  Azmy  PHI  problem  with  LD  and  an  S6  angular 
quadrature  are  shown  in  table  6.4.  The  spectral  radii  listed  in  the  table  were 
again  determined  using  the  linear  regression  of  the  slope  of  the  linearized  plots 


Table  6.4.  DI  with  LD  results. 


Cross  Section 

Mesh 

10 

20 

40 

80 

160 

10x10 

0.101 

0.115 

0.089 

0.059 

0.044 

20x20 

0.101 

0.117 

0.118 

0.075 

0.053 

40x40 

0.194 

0.099 

0.124 

0.094 

0.055 

80x80 

0.269 

0.216 

0.153 

0.140 

0.099 

160x160 

0.245 

0.302 

0.160 

0.179 

0.133 

Unlike  the  zeroth  spatial  moment  method,  the  convergence  rates  for  DI 
method  with  first  spatial  moment  methods  do  increase  slightly  with  larger 
problems.  However,  the  spectral  radius  is  still  much  better  than  the  DSA 
methods  and  the  DI  method  works  for  all  the  cross  sections  tested  (did  not 
diverge).  Also,  and  for  the  WDD  results  as  well,  the  calculated  spectral  radii  are 
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an  asymptotic  value  from  the  plots  of  convergence  rates.  Just  using  the  problem 
set  tolerance  of  10  6  would  have  given  lower  spectral  radii. 

The  DI  method  has  demonstrated  an  improved  performance  over  DSA  for 
this  particular  problem.  The  rate  of  convergence  for  zeroth  spatial  moments 
methods  is  clearly  superior  for  DI.  The  convergence  rate  stays  almost  constant 
while  the  DSA  method  increased  strongly  with  an  increase  in  the  number  of  cells. 
The  first  moment  methods  also  had  good  improvement  in  the  rate  of 
convergence,  and  the  DI  methods  were  able  to  solve  the  problem  for  cross 
sections  the  DSA  method  diverged  on. 

3.  Azmy  PHI  Timing  Analysis 

The  DI  method  showed  good  improvement  over  the  DSA  method 
performance  for  the  Azmy  PHI  problem,  particularly  for  the  zeroth  spatial 
moment  method  of  WDD.  Two  questions  to  be  answered  are:  where  does  the  DI 
method  spend  its  computational  effort;  and  how  does  the  effort  change  as  the 
number  of  cells  increase?  An  intrinsic  FORTRAN  timing  function  was  used  to 
determine  the  amount  of  time  spent  in  each  portion  of  the  DI  iteration. 

The  DI  iteration  is  separated  into  two  parts  for  timing  purposes;  discrete 
ordinates  sweep  cell  calculations;  and  the  partial  current  problem.  The  discrete 
ordinates  sweep  cell  calculations  will  be  further  broken  down  into  the  within  cell 
calculation  and  the  discrete  ordinates  sweep.  For  timing  purposes,  only  a  single 
within  cell  calculation  and  discrete  ordinates  sweep  will  be  timed.  The  actual 
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times  can  be  scaled  from  these  time  values.  The  partial  current  problem  is 
further  separated  into  the  collapsing  /  setting  up  the  partial  problem  and  the 
time  needed  for  the  CXML  library  routine  to  solve  the  partial  current  problem. 
The  timing  analysis  was  done  for  the  DI  method  with  WDD  and  LD  using  S6 
and  a  cross  section  O’  -10  . 

Zeroth  Spatial  Moment  Methods 

The  WDD  results  of  the  time  analysis  for  the  main  parts  of  a  DI  iteration: 
the  iterations  time;  among  cell  calculation  time;  and  partial  current  problem  time 
are  shown  in  table  6.5.  As  can  be  seen  in  the  table,  the  time  for  the  partial 
current  problem  is  most  of  the  iteration  time,  more  than  three  times  the  discrete 
ordinates  sweep  cell  calculation  time.  Additionally,  separate  log-log  plot  shows 
that  the  iteration  portions  of  the  code  scale  linearly  with  the  number  of  cells, 
with  a  slope  of  1.00.  SC  gave  nearly  identical  timing  results  for  the  zeroth 
spatial  moment  tests. 


Table  6.5.  WDD  Iteration  Timing. 


Time  (s) 

Number 

of  Cells 

Iteration 

Partial 

Current 

Problem 

Among  Cell 

Calculations 

100 

0.063 

0.047 

0.019 

400 

0.234 

0.188 

0.047 

1600 

0.938 

0.766 

0.172 

6400 

4.000 

3.281 

0.656 
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The  WDD  results  of  the  time  analysis  for  the  partial  current  problem  are 


shown  in  table  6.6.  As  can  be  seen  in  the  figure,  the  predominance  of  the  time 
for  the  partial  current  problem  is  used  for  collapsing  and  setting  up  the  partial 
current  problem  due  to  the  number  of  matrix  multiplications  that  are  done. 

Table  6.6.  WDD  Partial  current  problem  timing. 


Time  (s) 

Number 

of  Cells 

Partial 

Current 

Problem 

Collapsing 

Direct 

Solver 

100 

0.047 

0.031 

0.016 

400 

0.188 

0.141 

0.031 

1600 

0.766 

0.547 

0.172 

6400 

3.281 

2.234 

0.906 

The  WDD  results  of  the  time  analysis  for  the  discrete  ordinates  sweep  cell 
calculations  is  shown  in  table  6.7.  As  can  be  seen  in  the  table,  most  of  the  time 
is  for  within  cell  calculation,  again  doing  the  matrix  multiplications,  and  is  about 
twice  the  time  for  the  discrete  ordinates  sweep. 
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Table  6.7.  WDD  Discrete  ordinates  sweep  timing. 


Time  (s) 

Number 

of  Cells 

Among  Cell 

Calculations 

Update 

Scattering 

Source 

DO 

Sweep 

100 

0.019 

0.016 

0.004 

400 

0.047 

0.031 

0.016 

1600 

0.172 

0.109 

0.063 

6400 

0.656 

0.391 

0.266 

First  Spatial  Moment  Methods 

A  similar  analysis  was  done  for  the  LD  spatial  method.  The  LD  results  of 
the  time  analysis  for  the  main  parts  of  a  DI  iteration:  the  iterations  time; 
discrete  ordinates  sweep  cell  calculation  time;  and  partial  current  problem  time 
are  shown  in  table  6.8.  As  can  be  seen  in  the  table,  the  time  for  the  partial 
current  problem  is  most  of  the  iteration  time,  more  than  ten  times  the  among  cell 
calculations  time.  This  is  significantly  more  than  the  zeroth  spatial  moment 
methods,  and  due  to  the  additional  matrix  multiplications  used  in  collapsing  to 
set  up  the  partial  current  problem.  Additionally,  separate  log-log  plots  show  that 
the  iteration  portions  of  the  code  scale  linearly  with  the  number  of  cells  with  a 
slope  of  0.9969.  LC  gave  nearly  identical  timing  results  for  the  first  spatial 
moment  tests. 
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Table  6.8.  LD  Iteration  timing. 


Time  (s) 

Number 

of  Cells 

Iteration 

Partial 

Current 

Problem 

Among  Cell 

Calculations 

100 

0.469 

0.422 

0.047 

400 

1.828 

1.656 

0.172 

1600 

7.328 

6.688 

0.609 

6400 

29.547 

26.969 

2.484 

The  LD  results  of  the  time  analysis  for  the  partial  current  problem  time 
are  shown  in  table  6.9.  As  can  be  seen  in  the  table,  the  predominance  of  the 
time  for  the  partial  current  problem  is  collapsing  and  setting  up  the  partial 
current  problem.  The  CXML  direct  solver  actually  takes  the  same  amount  of 
time  as  the  zeroth  spatial  moment  methods.  This  is  to  be  expected,  the  actual 
problem  size  is  the  same  for  both  methods. 

Table  6.9.  LD  Partial  current  problem  timing. 


Time  (s) 

Number 

of  Cells 

Partial 

Current 

Problem 

Collapsing 

Direct- 

Solver 

100 

0.422 

0.422 

0.016 

400 

1.656 

1.625 

0.031 

1600 

6.688 

6.469 

0.172 

6400 

26.969 

25.875 

0.906 
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The  LD  results  of  the  time  analysis  for  the  discrete  ordinates  sweep  cell 


calculation  time  is  shown  in  table  6.10.  As  can  be  seen  in  the  table,  most  of  the 
time  is  for  the  within  cell  calculation  again,  updating  the  scattering  sources, 
significantly  more  than  the  discrete  ordinates  sweep. 

Table  6.10.  LD  Discrete  ordinates  sweep  timing. 


Time  (s) 

Number 

of  Cells 

Among  Cell 

Calculations 

Update 

Scattering 

Source 

DO 

Sweep 

100 

0.047 

0.031 

0.016 

400 

0.172 

0.141 

0.031 

1600 

0.609 

0.531 

0.078 

6400 

2.484 

2.156 

0.328 

The  timing  analysis  of  the  DI  method  showed  three  important  points. 
First,  the  problem  iteration  time  scales  linearly  with  the  number  of  cells.  This 
will  be  important  when  using  the  DI  method  to  solve  very  large  problems  in 
higher  dimensions.  Second,  the  among  cell  calculations  using  a  within  cell 
calculation  followed  by  a  discrete  ordinates  sweep  is  an  efficient  way  to  update 
the  cell  edge  values.  The  computation  cost  of  the  discrete  ordinates  sweep  cell 
algorithm  is  less  than  doing  two  within  cell  calculations  to  update  cell  edge 
values.  Lastly,  most  of  the  computational  effort  for  an  iteration  is  in  setting  up 
and  solving  the  partial  current  problem.  The  discrete  ordinates  sweep  cell 
calculations  are  a  smaller  part  of  the  computational  effort,  particularly  with  first 
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moment  methods.  The  timing  analysis  showed  that  additional  effort  could  be 
applied  to  improving  C,  without  a  significant  computational  cost  as  is  done  in 
the  discrete  ordinates  sweep  method. 

B.  Where  TSA  Fails 

The  Azmy  PHI  problem  showed  how  the  DSA  method  lost  effectiveness, 
or  converged  slowly  across  a  variety  of  cross  sections  and  meshes.  Another 
periodic  horizontal  interface  (PHI)  problem  reported  by  Chang  and  Adams  (5:  1) 
demonstrated  how  TSA  methods  diverged  for  certain  cross  section  combinations. 
This  next  problem  also  uses  pairs  of  cross  sections,  but  the  layout  and  source  are 
slightly  different  from  the  Azmy  PHI  used  in  the  last  section.  This  next 
problem,  hereafter  referred  to  as  the  Chang  problem(5:  11),  uses  a  fixed  mesh  of 
100  cm  by  200  cm  and  varies  the  two  different  cross  sections  during  these  tests. 
Each  cell  is  set  at  1  cm  by  1  cm  and  there  are  incident  boundary  currents  on  the 
bottom  and  left  sides  with  no  sources  within  the  problem.  The  layout  of  this 
problem  is  shown  in  figure  6.4. 
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aT  =  0 


Figure  6.4.  Problem  variables  for  the  Chang  Periodic  Horizontal  Interface 
(PHI). 

The  problem  was  done  for  different  scattering  ratios  and  to  a  tolerance  of 
10-7 using  the  cell  average  scalar  flux.  As  with  the  DSA  PHI  problem,  the 
measure  of  effectiveness  used  in  this  article  was  the  spectral  radius,  which  is 
determined  using  the  relation  described  in  equation  (6.1)  for  the  published  results. 
The  numerical  results  of  the  TSA  method  using  diamond  difference  (DD)  showed 
that  certain  cross  sections  caused  the  method  to  diverge,  as  can  be  seen  by 
spectral  radii  greater  than  one.  This  data  will  be  presented  later  in  this  chapter. 
The  DI  method  was  done  with  a  single  discrete  ordinates  sweeping  method 
initially  for  comparison.  The  DI  method  used  the  same  linear  regression 
procedure  that  was  used  with  the  Azmy  PHI  to  determine  the  convergence  rates. 
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1.  Weighted  Diamond  Difference  Comparison 


An  example  of  how  the  SRD  of  the  scalar  flux  changed  per  iteration  for  the  DI 


method  with  WDD  is  shown  in  figure  6.5  for  a  scattering  ratio  of  c  =  0.9. 
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Figure  6.5.  Convergence  rates  for  the  Chang  Periodic  Horizontal  Interface 
(PHI)  problem  with  DI  using  WDD  and  S6  for  a \  =10-4  and  <7-,  at 
various  cross  section  combinations  with  a  scattering  ratio  of  0.9. 
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As  can  be  seen  in  figure  6.5,  there  are  particular  cross  section 
combinations  for  which  the  problem  converges  slower.  The  fastest  convergence 
occurs  when  the  two  cross  sections  are  the  same,  making  it  a  homogeneous 
problem.  The  rates  of  convergence  are  listed  in  table  6.11  along  with  TSA  results 
(5:  11)  for  comparison. 
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Table  6.11.  Chang  PHI  Test  for  c  =  0.9,  with  DI 
WDD  using  S4  ,  and  TSA  DD  using  S4  /  S2  ■ 


cr1(=1.00 

(Tu= 10000.0 

Spectral 

Radius 

Spectral 

Radius 

Spectral 

Radius 

&2,t 

m 

TSA 

EH 

TSA 

EH 

TSA 

mm 

mm 

4.5255 

0.002214 

0.0397 

0.01 

0.212961 

O 

o 

SO 

OJ 

0.273905 

mm 

0.1745 

1.00 

0.273779 

mm 

0.4521 

100.0 

0.023046 

mm 

mm 

0.0783 

mm 

0.0398 

0.12314 

0.0089 

Note  that  for  four  different  cross  section  combinations,  the  TSA  method 
diverged.  This  is  indicated  by  a  spectral  radius  greater  than  one.  However,  for 
this  case,  the  DI  method  performed  well  having  a  spectral  radius  less  than  0.3  for 
all  cross  sections  and  only  using  a  single  discrete  ordinates  sweep  per  iteration. 
For  only  one  combination  of  cross  sections,  where  TSA  worked  well,  the  TSA 
spectral  radius  was  smaller  than  the  DI  method  spectral  radius.  Note  these 
problems  were  done  with  similar  spatial  methods  and  similar  angular 
quadratures. 

For  higher  scattering  ratios,  c  —  0.99  and  a  higher  order  angular 
quadrature,  the  following  comparisons  (5:  14)  can  be  made  in  table  6.12. 
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Table  6.12.  Chang  PHI  Test  for  c  =  0.99  with  DI  Results 
WDD  using  S6  ,  and  TSA  DD  Results  using  S6  /  S2- 


(7]  t  =  \  .00 

Spectral 

ladius 

,t 

DI 

TSA 

0.0001 

0.4858 

32.5264 

0.01 

0.4802 

20.1193 

1.00 

0.1017 

0.5865 

100.0 

0.0928 

1.2958 

10000.0 

0.1196 

0.5458 

Note  that  for  the  same  cross  section  combinations  that  diverged  in  the 
previous  problem,  the  TSA  method  also  diverges  for  this  case,  as  well  as  another 
combination  of  cross  sections.  The  DI  method  converges  for  this  problem  using 
one  discrete  ordinates  sweep  per  iteration.  The  performance  is  somewhat  slower 
for  the  particular  cross  section  pairs  where  TSA  diverged  for  this  case.  However, 
the  rates  of  convergence  for  the  other  cross  sections  remains  about  the  same  or 
slightly  faster  than  the  c  =  0.9  scattering  ratio  case  while  the  TSA  method  is 
much  slower. 

Although  TSA  method  did  not  give  results  for  scattering  ratios  of  c  =  1.0 , 
since  it  had  already  diverged  for  lower  scattering  ratios,  the  DI  method  was  also 
done  for  totally  scattering  problems  to  see  if  the  DI  method  would  solve  these 
problems  with  a  single  discrete  ordinates  sweep.  This  would  make  this  particular 
problem  as  difficult  as  it  could  be. 
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As  can  be  seen  in  table  6.13,  the  rate  of  convergence  for  this  scattering 


ratio  is  again  much  slower  for  certain  cross  sections  combinations. 

Table  6.13.  Chang  PHI  Test  for  c  =  1.0,  DI  with  WDD  using  S6 


,t 

<JU=  0.0001 

(J1(=1.00 

(Tj,  =  10000.0 

Spectral 

Radius 

Spectral 

Radius 

Spectral 

Radius 

1.00E-04 

0.006858 

0.653732 

0.420436 

1.00E-02 

0.234153 

0.641062 

0.406724 

1.00E+00 

0.651778 

0.111584 

0.04627 

1.00E+02 

0.397283 

0.038089 

3.7E-05 

1.00E+04 

0.395913 

0.038089 

2.12E-06 

For  the  combinations  where  TSA  diverged  previously,  the  DI  method 
performance  was  again  slower  but  still  converged.  The  other  cross  sections 
continued  to  converge  at  the  same  or  a  faster  rate. 

Again,  the  DI  results  with  WDD  are  presented  here  using  only  a  single 
discrete  ordinates  sweep  per  iteration  The  Chang  PHI  problem  was  challenging 
for  the  DI  method  for  certain  cross  section  combinations,  but  these  are  the  same 
cross  section  combinations  which  caused  the  TSA  method  to  diverge  for  this 
particular  problem. 

2.  Other  Spatial  Method  Comparison  for  Chang  PHI  Problem 

The  next  section  demonstrates  the  effect  the  adaptive  discrete  ordinates 
sweep  has  on  convergence  rates  for  this  problem.  The  adaptive  technique  is 
applied  to  the  Chang  PHI  problem  for  the  particular  cross  sections  that 
challenged  the  DI  method.  These  cross  sections,  <JU  =1.0  and  <7lt  =0.0001,  are 
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the  same  cross  sections  where  the  TSA  method  diverged  for  all  scattering  ratios. 


The  results  for  these  cross  sections  using  the  four  spatial  DI  methods  are  listed  in 
table  6.14  with  scattering  ratio  of  c  -  0.9  .  As  a  comparison  of  note,  the  TSA 
method  using  DD  had  a  spectral  radius  of  4.5255  or  diverged  (5:  11). 


Table  6.14.  Chang  PHI  Problem 
for  c=0.9  using  S6  . 


TSA  DD 

4.5255 

Spectral  Radius 

Method 

lx  DO 

Adaptive 

sweep 

DO  sweep 

WDD 

0.259 

0.198107 

SC 

0.431519 

0.25439 

LD 

0.466337 

0.249747 

LC 

0.473478 

0.243725 

All  the  spatial  methods  show  improvement  in  the  rate  of  convergence  for 
the  adaptive  discrete  ordinates  sweep  over  a  single  DO  sweep  per  iteration. 

The  same  case  was  done  again  for  a  scattering  ratio  of  c  =  0.99  .  The 
results  are  listed  in  table  6.15.  Again,  note  the  TSA  method  diverged  (5:  14)  for 
this  case.  The  adaptive  discrete  ordinates  sweep  technique  shows  improvement 
over  a  single  DO  sweep  per  iteration  in  the  rates  of  convergence  for  all  the  spatial 
methods  tested  with  the  DI  method. 
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Table  6.15.  Chang  PHI  Problem  for 
c=0.99  using  S6 . 


TSA  DD 

32.5264 

Spectral  Radius 

Method 

lx 

Adaptive 

DO  sweep 

DO  sweep 

WDD 

0.4858 

0.3907 

SC 

0.7217 

0.5471 

LD 

0.8285 

0.4851 

LC 

0.7783 

0.5058 

Although  the  TSA  method  was  not  done  for  a  scattering  ratio  of  c  =  1.0  , 

this  combination  of  cross  sections  caused  the  DI  method  with  SC  to  diverge  as 

well.  The  TSA  PHI  problem  was  done  again  using  the  adaptive  discrete 

ordinates  sweep  and  the  results  are  shown  in  table  6.16. 

Table  6.16.  Chang  PHI  Problem 
for  c=1.0  using  S6 


TSA 

N/A 

Spectral  Radius 

Method 

lx 

Adaptive 

WDD 

0.6537 

0.5786 

SC 

N/A 

0.9001 

LD 

0.9468 

0.7140 

LC 

0.8707 

0.6438 
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The  results  also  show  that  the  adaptive  discrete  ordinates  sweep  technique 
improve  the  rates  of  convergence  for  the  DI  method.  The  adaptive  DO  sweep 
technique  also  stabilizes  the  SC  spatial  method,  which  had  previously  diverged 
for  this  problem. 

This  section  showed  that  the  improved  performance  for  the  DI  method 
over  TSA  for  this  particular  problem.  The  DI  WDD  method  converged  reliably 
where  TSA  DD  did  not.  The  DI  method  rate  of  convergence  was  considerably 
faster  than  TSA  when  TSA  did  work. 
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VII.  Additional  Tests 


This  section  shows  the  development  and  analysis  of  the  adaptive  DO 
sweep  technique  that  was  used  in  chapter  six.  Earlier  testing  on  the  Chang  PHI 
problem  showed  areas  where  the  DI  method  performance  was  challenged  for 
certain  spatial  methods  and  scattering  ratios.  In  an  attempt  to  fully  challange 
the  method,  another  degree  of  interfaces  or  a  checkerboard  of  alternating  cells 
was  added.  Also,  to  further  stress  the  DI  method  with  this  problem,  the  cross 
sections  that  caused  the  TSA  method  to  diverged  and  showed  slower  convergence 
rates  for  DI  were  chosen.  The  cross  section  values  used  were  <7,  t  =1.0  and 
(J2l  =  0.0001 ,  and  the  scattering  ratio  is  set  to  one. 

A.  The  Checkerboard  Problem 

The  problem  used  incident  currents  on  the  left  and  bottom  side,  like  the 
Chang  PHI  problem.  A  diagram  of  the  problem  is  shown  in  figure  7.1.  For  this 
problem,  the  S6  angular  quadrature  was  used.  The  number  of  x  and  y  cells 
were  the  same  for  each  case  and  varied  from  25  to  125  each. 
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Figure  7.1.  Problem  variables  for  the  checkerboard  problem. 

1.  Spatial  Method  Performance 

As  was  shown  in  the  chapter  six,  other  spatial  methods  did  not  perform  as 
well  as  WDD  for  the  Chang  PHI  problem  with  these  particular  cross  sections  and 
scattering  ratio.  The  rate  of  convergence  for  the  checkerboard  problem  can  be 
seen  in  figure  7.2  for  all  the  spatial  methods  currently  implemented  in  DI  using  a 
single  discrete  ordinates  sweep. 
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Figure  7.2.  DI  Convergence  rates  for  different  spatial  methods  on  the 
checkerboard  problem  for  a  single  DO  sweep. 


As  can  be  seen  in  figure  7.2,  SC  does  not  converge  for  this  problem  (like 
the  Chang  PHI  problem)  with  one  DO  sweep.  The  first  spatial  moment  methods 
are  slow  to  converge  for  tolerance  >  10-4  while  WDD  converges  in  the  fewest 
number  of  iterations.  The  performance  of  the  spatial  methods  for  the  Chang  PHI 
problem  gave  similar  results.  This  performance  in  figure  7.2  for  the  DI  method 
indicated  that  the  checkerboard  problems  were  taxing  the  DI  method  using  only 
one  DO  sweep  per  iteration  for  several  spatial  methods. 

2.  Edge  Distribution  Improvement 

The  initial  attempt  at  how  the  number  of  DO  sweeps  per  iteration 
influences  the  overall  performance  is  shown  in  figure  7.3.  In  this  case,  the 
number  of  discrete  ordinate  sweeps  was  increased  from  one  time  per  iteration  to 
three  times  per  iteration  for  each  of  the  four  spatial  methods. 
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Figure  7.3.  DI  Convergence  rates  for  different  spatial  methods  on  the 
checkerboard  problem  for  a  three  DO  sweeps. 


Figure  7.3  shows  several  important  points.  With  the  additional  discrete 
ordinate  sweeps,  SC  is  now  converging  slowly  as  opposed  to  diverging,  and  the 
first  spatial  moment  methods  are  converging  faster,  especially  LD.  However,  the 
effect  on  WDD  is  small,  there  is  little  change  in  the  convergence  rate. 

This  spatial  method  dependence  led  to  the  concept  of  letting  the  code 
decide  how  much  effort  to  put  into  discrete  ordinate  sweeps,  or  an  adaptive 
technique  to  estimate  how  many  discrete  ordinates  sweeping  cell  calculations  to 
do.  For  spatial  methods  that  are  working  well,  like  WDD,  there  is  little 
advantage  to  doing  additional  discrete  ordinate  sweeps.  For  spatial  methods  that 
are  not  performing  well,  like  SC,  more  effort  in  the  discrete  ordinates  sweeping 
cell  calculations  should  help  the  problem  converge  quicker.  There  should  be  a 
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limit  to  the  maximum  number  of  discrete  ordinates  sweeping  cell  calculations.  A 
first  estimate  of  what  the  maximum  should  be  is  based  on  the  results  of  the 
timing  analysis.  The  partial  current  problem  time  during  an  iteration,  for  the 
first  spatial  moment  methods,  was  about  ten  times  longer  than  the  discrete 
ordinates  sweeping  cell  calculations.  Thus  ten  was  chosen  as  the  maximum,  as  it 
would  at  most  double  the  iteration  time  for  the  first  spatial  moment  methods. 

The  adaptive  technique  used  the  ratio  of  the  maximum  value  in  the  SRD  of  the 
edge  distribution  f  for  the  current  and  previous  iteration. 

..  Max(SRD(C‘ 

Number  of  DO  sweeps  =  Maxvaluex -  —  (7.1) 

Max(SRD(£  ,Q  )) 

The  ratio  of  the  maximum  values  of  the  SRD  of  f  should  be  less  than  one  for  a 
method  that  is  converging  and  much  less  for  one  that  is  converging  quickly.  For 
methods  that  are  working  well,  only  one  discrete  ordinates  sweeping  calculation  is 
enough  to  improve  the  estimate  of  tf  and  the  ratio  should  reflect  that.  Methods 
that  need  additional  effort  would  have  more  discrete  ordinate  sweeps  up  to  ten. 
The  results  of  the  adaptive  algorithm  are  shown  in  figure  7.4. 
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Figure  7.4.  DI  Convergence  rates  for  different  spatial  methods  on  the 
checkerboard  problem  for  the  adaptive  DO  sweep  technique. 

As  in  the  multiple  calculations  shown  in  figure  7.3,  the  adaptive  technique 
used  in  figure  7.4  shows  similar  results  with  a  few  key  differences.  The 
performance  of  the  first  spatial  moment  methods,  LD  and  LC,  is  better, 
converging  in  fewer  iterations.  Also  performance  of  SC  which  was  slow  or  even 
diverged  in  the  previous  two  cases,  is  much  improved.  The  SC  method  now 
converges  readily.  Again,  the  improvement  of  WDD  is  not  significant,  it  has 
been  working  well  previously.  The  number  of  iterations  need  to  reach  a  tolerance 
of  1CT7  is  shown  in  table  7.1. 
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Table  7.1.  Number  of  DI  iterations  for  the  checkerboard  problem. 


DO  Sweep  Cell  Calculation  Technique 

Spatial 

Method 

lx 

3x 

Adaptive 

WDD 

28 

23 

21 

SC 

Div 

>100 

57 

LD 

>100 

55 

35 

LC 

54 

38 

29 

The  improvement  in  performance  for  some  spatial  methods  is  considerable 
for  this  simple  adaptive  technique.  Another  optimization  of  the  technique,  or 
different  choice  for  the  maximum  number,  may  give  even  better  performance. 
However,  this  simple  adaptive  method  is  sufficient  to  show  the  robustness  of  the 
DI  method  and  the  importance  of  efficiently  improving  the  estimate  of  C,  values 
for  difficult  problems. 

Table  7.2.  Total  number  of  DO  sweeps  for  the  checkerboard  problem 


Number  of  DO  Sweep  Calculations 

Spatial 

Method 

lx 

3x 

Adaptive 

WDD 

28 

69 

83 

SC 

Div 

>300 

393 

LD 

>100 

165 

220 

LC 

54 

114 

156 

Table  7.2  shows  the  total  number  of  discrete  ordinate  sweeps  done  for 
each  of  the  different  spatial  methods.  For  each  case,  the  total  number  of  discrete 
ordinate  sweeps  is  two  to  three  times  more.  However,  these  calculations  are  the 
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comparatively  inexpensive  part  of  the  overall  iteration  calculation.  Timing 
analysis  shows  that  setting  up  and  solving  the  partial  current  problem  is  most  of 
the  computational  time  of  an  iteration.  Improving  the  estimate  of  £  values 
during  the  among  cell  calculations  is  more  important  than  setting  up  and  solving 
more  partial  current  problems. 

3.  Timing  Analysis 

The  maximum  number  of  discrete  ordinate  sweeps  was  chosen  so  as  to  at 
most  double  the  iteration  time  for  the  first  spatial  moment  methods.  This  leads 
to  the  question  of  what  does  the  additional  calculations  do  to  the  total  time  to 
solving  the  problem?  The  plot  of  the  total  time  for  the  DI  WDD  method 
comparing  the  difference  in  time  as  the  number  of  cells  increase  is  shown  in  figure 
7.5. 

As  can  be  seen  in  figure  7.5,  the  total  time  is  slightly  more,  even  though 
the  number  of  discrete  ordinate  sweeps  cell  calculations  tripled  as  shown  in  table 
7.2.  This  is  due  to  the  fewer  number  of  total  number  of  iterations  shown  in  table 
7.1  which  offset  the  time  for  the  additional  discrete  ordinate  sweeps  cell 
calculations. 
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Figure  7.5.  Comparison  of  total  convergence  time  for  the  checkerboard 
problem  with  DI  using  WDD  for  a  single  DO  sweep  and  the  adaptive  DO 
sweep  technique. 

The  first  spatial  moment  time  is  more  interesting,  as  the  maximum 
number  of  sweeps  was  chosen  for  these  methods  in  particular.  The  plot  of  the 
total  time  for  the  DI  LD  method  comparing  the  difference  in  time  as  the  number 
of  cells  increase  is  shown  in  figure  7.6. 
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Figure  7.6.  Comparison  of  total  convergence  time  for  the  checkerboard 
problem  with  DI  using  LD  for  a  single  DO  sweep  and  the  adaptive  DO 
sweep  technique. 

As  can  be  seen  in  the  figure,  the  total  time  is  consistently  less  for  the 
adaptive  technique  by  as  much  as  a  factor  of  two,  even  though  the  number  of 
discrete  ordinate  sweeps  cell  calculations  doubled  as  shown  in  table  7.2.  Again, 
the  decrease  in  time  is  due  to  the  fewer  number  of  total  iterations,  as  shown  in 
table  7.1,  which  offset  the  time  for  the  additional  discrete  ordinate  sweeps  cell 
calculations. 

The  checkerboard  problem  demonstrated  the  robustness  of  the  DI  method. 
The  problem  challenged  the  method  initially,  as  some  of  the  PHI  problems  did 
for  synthetic  acceleration  methods.  The  checkerboard  problem  however,  showed 
a  way  to  both  stabilize  diverging  spatial  methods  and  improve  the  convergence 


•  lx  Among 
O  Adaptive  Among 


1000 


10000 


Number  of  Cells 


148 


rates  methods  that  were  converging  slowly.  Furthermore,  the  new  technique  does 
not  come  at  a  computational  cost  penalty,  it  actually  improves  the  overall  speed 
of  the  method  for  first  spatial  moment  methods. 

B.  Scattering  Ratio  Horizontal  Interface  Problem 

Another  problem  that  was  tested  is  also  a  periodic  horizontal  interface, 
but  where  the  scattering  ratio  is  varied  between  layers  as  opposed  to  the  total 
cross  sections.  The  spatial  mesh  is  40  cells  by  40  cells  with  Ax  =  Ax  =  1.0.  A 
description  of  the  problem  is  shown  in  figure  7.13  for  a  total  cross  section  <7  = 1 . 


aR  =  0 

<7  =  1.0 

□  Region  A 

□  Region  B 
[~~~1  Source  Region 


Figure  7.7.  Problem  variables  for  the  scattering  ratio  PHI  problem. 


Each  region’s  scattering  ratio  are  systematically  changed  and  the  rates  of 
convergence  are  checked  for  total  cross  sections  of  0.1  ,  1.0,  and  10.0.  The 
problem  is  intended  to  stress  the  DI  method  for  the  first  moment  methods  by 
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creating  regions  where  the  current  along  an  ordinate  is  not  continuous  and  may 


produce  negative  current  artifacts. 

Spatial  Methods 

The  DI  results  of  the  scattering  ratio  PHI  for  the  total  cross  section  of  0.1 
is  shown  in  figure  7.8.  All  spatial  methods  converged  readily  for  the  scattering 
ratio  of  1.0  and  0.0  for  regions  A  and  B  respectively.  The  various  combinations 
of  scattering  ratios  are  shown  later  for  WDD  at  this  total  cross  section  in  table 
7.3. 
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Figure  7.8.  Convergence  rates  for  different  spatial  methods  versus 
iterations  using  S6  to  the  scattering  ratio  PHI  problem.  Total  cross 
section  is  0.1  and  scattering  ratios  of  1.0  and  0.0  for  regions  A  and  B 
respectively. 
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The  results  of  the  scattering  ratio  PHI  for  the  total  cross  section  of  1.0  is 


shown  in  figure  7.9.  Again,  all  spatial  methods  converged  readily  for  the 
scattering  ratio  of  1.0  and  0.0  for  regions  A  and  B  respectively.  The  various 
combinations  of  scattering  ratios  are  shown  later  for  WDD  at  this  total  cross 
section  in  table  7.4. 
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Figure  7.9.  Convergence  rates  for  different  spatial  methods  versus 
iterations  using  S6  to  the  scattering  ratio  PHI  problem.  Total  cross 
section  is  1.0  and  scattering  ratios  of  1.0  and  0.0  for  regions  A  and  B 
respectively. 


The  results  of  the  scattering  ratio  PHI  for  the  total  cross  section  of  10.0  is 
shown  in  figure  7.10.  As  with  the  previous  cases,  all  spatial  methods  converged 
readily  for  the  scattering  ratio  of  1.0  and  0.0  for  regions  A  and  B  respectively. 
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The  various  combinations  of  scattering  ratios  are  shown  later  for  WDD  at  this 


total  cross  section  in  table  7.4. 
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Figure  7.10.  Convergence  rates  for  different  spatial  methods  versus 
iterations  using  S6  to  the  scattering  ratio  PHI  problem.  Total  cross 
section  is  10.0  and  scattering  ratios  of  1.0  and  0.0  for  regions  A  and  B 
respectively. 


1.  Weighted  Diamond  Difference  Performance 

The  various  combinations  of  scattering  ratios  for  regions  A  and  B  are 
shown  in  table  7.3  for  a  total  cross  section  of  0.1.  For  combinations  where  both 
regions  that  are  totally  absorbing  or  with  a  scattering  ratio  of  0.0,  the  DI  method 
converged  in  one  iteration.  As  with  the  previous  convergence  rates,  the  values 
are  determined  by  a  linear  regression  of  the  linearized  maximum  SRD  plots 
described  in  chapter  six.  In  general,  the  convergence  rates  increase  slightly  for 
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higher  scattering  ratios,  but  the  overall  performance  is  good.  For  the  other 

spatial  methods,  the  convergence  rates  were  similar. 

Table  7.3.  WDD  results  for  scattering  ratio  PHI  problem  with  total  cross 
section  ex  =  0. 1  as  scattering  ratio  varies. 


c  Region  A 

c  Region  B 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.0 

0 

0.0822 

0.0772 

0.0630 

0.1501 

0.1429 

0.2 

0.00002 

0.0518 

0.1104 

0.0622 

0.1338 

0.1385 

0.4 

0.0656 

0.0759 

0.0624 

0.1225 

0.1666 

0.1574 

0.6 

0.1051 

0.0629 

0.1464 

0.1570 

0.1831 

0.1957 

0.8 

0.0780 

0.1335 

0.1588 

0.1660 

0.1864 

0.1892 

1.0 

0.1119 

0.1585 

0.1529 

0.1449 

0.1992 

0.2953 

The  various  combinations  of  scattering  ratios  for  regions  A  and  B  are 
shown  in  table  7.4  for  a  total  cross  section  of  1.0  and  the  DI  method  with  WDD. 
Again,  where  both  regions  that  are  totally  absorbing  or  with  a  scattering  ratio  of 
0.0,  the  DI  method  converged  in  one  iteration.  Again,  the  convergence  rates 
increase  slightly  for  higher  scattering  ratios,  but  the  overall  performance  is  very 
good.  The  convergence  rates  are  slightly  better  that  the  previous  total  cross 
section  of  0.1  shown  in  table  7.3.  Again,  the  other  spatial  methods  had  similar 
performance  for  this  total  cross  section,  the  LD  method  is  presented  next. 
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Table  7.4.  WDD  results  for  scattering  ratio  PHI  problem  with  total  cross 
section  (7  =  1.0. 


c  Region  A 

c  Region  B 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.0 

0 

0.0161 

0.0643 

0.0681 

0.0685 

0.0657 

0.2 

0.00004 

0.0203 

0.0674 

0.0643 

0.0476 

0.0659 

0.4 

0.0499 

0.0506 

0.0532 

0.0403 

0.0439 

0.0640 

0.6 

0.0344 

0.0411 

0.0348 

0.0341 

0.0457 

0.0686 

0.8 

0.0280 

0.0315 

0.0413 

0.0527 

0.0542 

0.0875 

1.0 

0.0499 

0.0506 

0.0595 

0.0657 

0.0587 

0.1046 

2.  Linear  Discontinuous  Performance 

The  various  combinations  of  scattering  ratios  for  regions  A  and  B  are  shown  in 
table  7.5  for  a  total  cross  section  of  1.0  and  the  DI  method  with  LD.  Again, 
where  both  regions  that  are  totally  absorbing,  the  DI  method  converged  in  one 
iteration. 


Table  7.5.  LD  results  for  scattering  ratio  PHI  problem  with  total  cross 
section  o  -  1.0  .  Shading  indicates  strictly  positive  solutions. 


c  Region  A 

c  Region  B 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.0 

0 

2.0E-06 

0.0943 

0.0049 

0.0185 

0.1489 

0.2 

9.8E-8 

4.0E-06 

0.0005 

0.0858 

0.0685 

0.1518 

0.4 

5.1  E-7 

0.0003 

0.0680 

0.0906 

0.0823 

0.1548 

0.6 

0.0004 

0.0772 

0.1058 

0.0691 

0.0968 

0.1332 

0.8 

0.1156 

0.1009 

0.0738 

0.1005 

0.1114 

0.1773 

1.0 

0.0923 

0.1057 

0.1022 

0.1427 

0.1825 

0.2571 

154 


The  first  moment  methods,  LD  and  LC,  are  not  positive  methods  and  can 


return  negative  values  for  certain  problem  values.  A  test  was  done  to  see  if 
either  the  cell  scalar  flux  or  edge  currents  (and  hence  distributions)  were  negative 
during  any  iteration  for  the  solutions  presented  in  table  7.5.  The  results  are 
shown  in  table  7.5,  and  the  positive  values  are  indicated  by  the  shaded  cells.  The 
table  shows  that  the  LD  method  did  indeed  return  negative  values  for  most  of 
the  cases,  where  the  scattering  sources  less  than  0.8.  However,  for  all  the 
combinations  in  cross  sections,  the  method  was  still  able  to  converge  in  spite  of 
the  spatial  method  negative  artifacts.  The  DI  method  is  able  to  tolerate  some 
negative  values  as  demonstrated  by  this  case. 

The  various  combinations  of  scattering  ratios  for  regions  A  and  B  from 
figure  7.10  are  shown  in  table  7.6  for  a  total  cross  section  of  10.0  and  the  DI 
method  with  WDD.  In  general,  the  convergence  rates  are  constant  across  the 
range  of  scattering  ratios.  The  overall  convergence  rates  are  fast  for  WDD  and 
SC  had  similar  performance.  However,  for  this  total  cross  section  both  the  LD 
and  LC  spatial  method  did  not  converge  for  certain  scattering  ratios.  The 
negative  artifacts  returned  by  these  spatial  methods  prevented  the  DI  method 
from  working. 
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Table  7.6.  WDD  results  for  scattering  ratio  PHI  problem  with  total  cross 
section  (7  =  10.0. 


c  Region  A 

c  Region  B 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.0 

0 

5.9E-06 

0.0067 

0.0115 

0.0079 

0.0082 

0.2 

6.5E-6 

2.3E-06 

0.0063 

0.0135 

0.0143 

0.0084 

0.4 

0.0070 

0.0070 

0.0195 

0.0167 

0.0151 

0.0085 

0.6 

0.0126 

0.0157 

0.0162 

0.0162 

0.0153 

0.0087 

0.8 

0.0139 

0.0149 

0.0148 

0.0154 

0.0168 

0.0091 

1.0 

0.0077 

0.0083 

0.0080 

0.0082 

0.0085 

0.0067 

The  scattering  ratio  periodic  horizontal  interface  problem  also 
demonstrated  the  robustness  of  the  positive  spatial  methods  for  problems  that 
have  a  difference  in  scattering  ratio  at  cell  boundaries.  The  positive  methods 
performed  well  across  the  entire  range  of  cross  section  combinations  and  different 
of  total  cross  sections.  This  problem  also  highlighted  the  issue  of  how  the  DI 
method  responds  to  the  negative  artifacts  created  by  the  first  moment  methods. 
The  DI  method  is  able  to  tolerate  some  negative  values,  but  other  cases  will 
cause  it  to  fail.  An  obvious  approach,  is  to  refine  the  mesh  in  an  attempt  to  keep 
the  first  moment  methods  positive.  Another  commonly  used  approach  is  to 
impose  a  fix-up  and  set  the  negative  values  to  zero.  Both  these  approaches 
create  issues  for  the  rate  of  convergence.  However,  rather  than  address  this  issue, 
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it  is  more  desirable  to  devise  spatial  methods  that  are  strictly  positive  and  not 
contain  non-physical  numerical  artifacts. 
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VIII.  Conclusions  and  Recommendations 


This  research  showed  that  the  distribution  iteration  method  is  a  practical 
alternative  to  current  methods  and  suggests  possibilities  for  new  approaches  to 
solving  the  discrete  ordinates  system  of  equations. 

A.  Conclusions 

My  objectives  have  been  achieved.  The  distribution  iteration  method  was 
extended  to  2-d  Cartesian  Geometry  (objective  1)  and  demonstrated  using 
multiple  spatial  and  angular  quadratures,  including  quadratures  that  correctly 
meet  diffusion  limits  (objective  2).  Unlike  the  synthetic  acceleration  methods,  a 
different  derivation  is  not  required  to  change  spatial  methods,  as  this 
demonstrated.  The  global  problem  was  recast  as  a  finite- volume  particle 
conservation  formulation  (objective  3)  by  using  partial  currents,  rather  than 
partial-range  angular  integrals  of  the  directional  flux,  creating  the  global  partial 
current  problem.  This  change  was  not  only  shown  to  converge  in  fewer 
iterations,  but  also  provides  a  clear  methodology  for  the  extension  to  higher 
dimensions.  The  global  problem  is  defined  in  terms  of  only  the  spatial  average  of 
the  two  partial  currents  through  each  cell  face.  (Alternative  schemes  could 
include  higher  spatial  moments  and/or  cell  spatial  moments  as  well,  but  this 
would  increase  the  size  of  the  global  problem.)  Thus,  the  method  minimizes  the 
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size  of  the  global  problem  when  applied  to  higher-order  linear  spatial  quadratures 
(objective  4). 

The  distribution  iteration  method  efficiently  solved  problems  where  the 
synthetic  acceleration  methods  either  failed  or  lost  effectiveness;  and  testing  on 
more  challenging  problems  also  demonstrated  the  success  of  the  distribution 
iteration  method  (objective  5).  Despite  the  comparisons  with  synthetic 
acceleration  methods  in  this  research,  the  distribution  iteration  method  is  not  an 
acceleration  technique. 

PARDISO  was  evaluated  (objective  6)  and  found  to  be  extremely  efficient 
throughout  my  testing.  It  required  only  a  small  fraction  of  the  run  time  of  the 
code.  The  red/black  scheme  maximizes  opportunity  for  parallelization  (objective 
7a);  while  the  sweep  scheme  enhances  serial  performance  by  requiring  fewer  cell 
calculations  (objective  7b).  The  desirable  properties  of  the  method  that  constitute 
the  goals  for  the  research  have  been  nearly  fully  achieved  (objective  8): 

Robustness  -  the  method  has  been  demonstrated  for  a  broad  range  of  cross 
sections  and  scattering  ratios.  Convergence  was  sometimes  slow  or  divergent  for 
some  spatial  quadratures  but  an  inner  loop  with  an  adaptive  number  of  sweeps 
per  global  solver  call  (1  to  10)  offset  this  limitation.  This  is  the  one  area  that 
needs  future  work:  finding  a  better  sweep  scheme  for  updating  the  angular 
distributions  at  cell  faces. 
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Flexibility  -  Several  spatial  and  angular  quadratures  were  used.  These 


changed  the  numbers  in  the  matrices  for  the  cell  and  global  problems,  but  did 
not  require  changes  to  the  algorithm  (for  given  spatial  moments  carried). 

Parallelizability  -  The  global  problem  is  solved  by  the  PARDISO  solver, 
which  is  commercially  available  for  many  parallel  computing  systems.  The  sweep 
method  of  improving  angular  distributions  at  cell  faces  has  limited 
parallelizability,  but  the  red/black  alternative  also  demonstrated  here  is  ideal  for 
parallel  computation. 

Extensibility  -  the  method  could  readily  be  extended  to  3d.  The  algorithm 
design  is  unaffected,  only  changes  in  implementation  (array  dimensions,  cell 
indexing  and  translation  to  sparse  array  data  structure,  etc.)  are  needed. 

B.  Recommendations 

The  testing  showed  that  a  large  portion  of  the  code  run  time  was  for  the 
matrix  multiplications  required  in  different  steps.  Additional  time  savings  are 
likely  by  using  optimized  (vectorized)  matrix  multiplication  routines  and  a  more 
compact  data  structure  within  the  code. 

The  sweeping  scheme  needs  to  converge  faster  for  some  challenging 
problems,  such  as  the  checkerboard  problem.  This  might  be  achieved  by  tuning 
the  number  of  sweeps  per  global  solution,  by  applying  a  convergence  accelerator 
to  this  inner  iteration,  or  by  trying  variants  on  exactly  how  the  sweeps  are  done. 
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The  distribution  iteration  is  not  an  accelerator  for  the  source  iteration.  It 


may  be  improved  by  applying  some  acceleration  scheme(s);  this  is  an  open 
question. 
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Appendix  A:  Linear  Discontinuous  Equations  in  Slab  Geometry 

Equation  Section  1 

In  this  section,  the  relations  for  the  current  representation  of  the  linear 
discontinuous  method  in  slab  geometry  are  developed.  As  with  the  zeroth  spatial 
moment  methods,  the  usual  representation  found  is  in  terms  of  the  angular  flux. 
The  equations  for  the  outgoing  angular  flux  i//nghl  in  terms  of  the  incoming 
angular  flux  y/left ,  scattering  within  the  cell  SA ,  SA  and  emissions  EA ,  are 
previously  derived  and  presented  in  the  literature  (8:  222-223).  The  linear 
discontinuous  relation  for  the  edge  value  is: 


YRight  =yA +YX , 


(A.l) 


This  can  be  substituted  into  the  zeroth  and  x  moment  cell  balance 


equations: 


y/Right  -  y/Lefi  +ey/A=SA^i 

V 

3(y/Righl  +  y/Left  -  2 y/A)  +  ey/x  =SX 

n 


(A.2) 

(A.3) 


where  the  following  relation  is  defined: 


G  —  6  +  Ae  +  £ 


(A.4) 


and 


£  =  - 


crAr 

H 


(A.5) 


The  desired  relations  in  the  angular  flux  relation  are  found  using  equations  (A.l) 
through  (A.4): 
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,  Right  _6-2e  Left  Ax(6  +  e)  A  ,  A xe  x  ,  Ax(6  +  e)  A 

If/  — - If/  ~\  O  “t  O  \  I - , - Hi  , 

a  a  jU  a  jU  a\ju\ 


_  6  +  £  Left  |  Ax(3  +  g)^  Ax  Sx  I  A*(3  +  g) 
a  a  jU  a  fl  a|//| 


and 


W*  +^SA  +  +  3A* 

a  a\jA\  a  fl  a|//| 


(A.6) 

(A.7) 

(A.8) 


The  definition  for  the  current  jRigln  =  \[A\\f/Right  and  similarly  for  the  left,  allows 
the  transition  from  an  angular  flux  to  a  current  representation.  The  equations 


for  the  outgoing  quantities  are: 

■Right  _  6-2e  Left  |  A x(6  +  s)  s A  +^sx  l  A x(6  +  s)  £a 


a  _  6  +  g  j Left  +  Ax(3  +  g)  ^a  _  Ax  +  Ax(3  +  g)  pA 
a  jU  a  /il  a  jU  a  \ 


and 


x  _  ~3g  .Left  t  3A x  ^a  t  Ar(l  +  g)^,x  |  3A x 


V"  =~r\J 

a\M\ 


ju\  a \ju 

These  are  the  relationships  used  in  chapter  3. 


Ax 

!W 


(A.9) 

(A.10) 

(A.ll) 
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Appendix  B:  Zeroth  Spatial  Moment  Methods  Current  Equations  in 

XY-Geometry 

Equation  Section  2 

Step  Characteristic  Equations 

In  this  section,  the  relations  for  the  current  representation  of  the  step 
characteristic  method  are  developed.  The  usual  representation  for  the  cell 
relations  are  in  terms  of  the  angular  flux.  In  the  scaled  rectangular  cell  as  shown 
in  figure  A.l,  the  equations  for  the  outgoing  angular  fluxes  y/'op  and  l/f"s'u  in 
terms  of  the  incoming  angular  fluxes  y/bo"om  and  y/left  ,  scattering  within  the  cell 
SA  ,  and  emissions  E A ,  are  previously  derived  and  presented  in  the 
literature  (12:  21). 


Bottom 

Figure  B.l.  Rectangular  cell  for  zeroth  spatial  moment  methods.  Cell 
shows  problem  variables  used  for  the  discrete  ordinates  equations. 


The  angular  flux  relations  for  the  cell  shown  in  figure  B.l  in  the  step 
characteristic  spatial  quadrature  are: 
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WTop  =  aM0(£y)iisLeft  +  (l-a)e-£>rBottom  + 

— [(1  -  a)M0  (e  )  +  aM ,  (e  )\EA , 

Tj 


(B.l) 


y/Right  =  MJ£v)YBottom +—MAev)SA  +^tMAev)Ea  , 

y  jj  )  \rl\ 


Ay 

n 


(B.2) 


Va  =  aM0(£y)i//Left  +  [(1  -  a)M0(£y)  +  aMl(£y)]y/DO“om  + 

^{\-a)Ml{£y)  +  aM2(£y)]SA+^[{\-a)Ml{£y)  +  aM2{£}y\EA.  ( 

Here  Ll  and  Tj  are  the  direction  cosines  along  the  x  and  y  axis  respectively  from 


Bottom 


the  angular  quadrature,  £  =  is  the  optical  thickness  in  the  y  direction, 

y  11 


a  =  —  is  a  parameter  for  the  cell,  and  M0(£  )  ,  M^(£  ),  and  M2(£v)  are  the 


exponential  moment  functions  (9:  27). 

The  definition  for  the  currents,  jRlght  =| ju\y/Rlght  and  jTop  =  \r/\i//Top  and 
similarly  for  the  left  and  bottom,  allows  the  transition  from  angular  fluxes  to  a 
current  representation.  The  equations  (B.l)  through  (B.3)  in  a  current 
representation  are: 


jTop  MaM0(£y)jLefi  +(l-a)e-^jBottom  +Ay[(\-a)MQ(£y)  +  aM](£y)]SA  + 
Av[(  1  -  a)M0  {£  )  +  aMx{£  )]EA, 


(B.4) 


■Right  _  J/([ 

N 


J'  °  =T~jM0(£  )jB°tt0m + 


/'I  A v 

v 


M\(£y)S  + 


\M  Ay 

V 


M\(£y)E  , 


(B.5) 
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¥a  =  ^aM{){£y)jLel1  +  ^|[(1  -  a)M{){ey)  +  aMx  (£y)\jBollom  + 

—[{l-a)Ml(£y)  +  aM2(ey)]SA+—[{\-a)Ml(£y)  +  aM2(ey)\EA. 

rj  *  s  jj  >  > 


(B-6) 


These  are  the  relations  presented  in  chapter  four. 


Weighted  Diamond  Difference  Equations 


In  the  rectangular  cell  as  shown  in  figure  B.l,  the  equations  are  again 
found  for  the  outgoing  angular  fluxes,  y/,op  and  y/"ght ,  in  terms  of  the  incoming 
angular  fluxes  y/ho“om  and  y/kft ,  scattering  within  the  cell  SA ,  and  emissions  EA  . 
The  WDD  relations  begin  with  the  cell  balance  equation  (3:  215): 


-  —  —  y/^efl  )  +  H—(ij/ToP  _  yjBottom ^  _j_  (jy/A  =  gA  +  £A, 

Ax  Ay 


(B.7) 


with  the  weighted  diamond  difference  assumption 


A  =  1  +  0 ^  )  Right  J  l_0^_  Left 
2  2 


(B.8) 


l  +  ay 


1  -  ay 


Bottom 


(B.9) 


The  weights,  ax  and  ay ,  change  the  relations  from  a  diamond  difference 
to  step  spatial  quadrature  using  the  following  relations: 


ax  =  coth(— )-  — , 
2  e’ 


(B.10) 


ay  -  coth(— . 
2 


(B.ll) 


166 


Again,  fi  and  //  are  the  direction  cosines  along  the  x  and  y  axis  respectively  from 


the  angular  quadrature,  £y  =~^~  is  the  optical  thickness  in  the  y  direction, 


\n\ 


and  £ 


<jAx 


v  .  .  is  the  optical  thickness  in  the  x  direction.  These  relations  are 

H 


numerically  ill  conditioned  for  optically  thin  or  thick  cells,  but  may  be 
equivalently  expressed  using  exponential  moment  functions  (9:  27)  as 


«'  =  2p(ex\ 


(B.12) 


where: 


P(ex)  ■ 


M\(£x) 


Mq(£x) 

and  similarly  for  the  y  component.  A  new  notation  for  the  weights  can  be 
written  as: 


(B.13) 


w  _  1-2 p(£x) 


ox  _  1  +  2  p{sx) 
°Out  ^ 


(B.14) 

(B.15) 


The  definition  for  the  currents  jR‘gh>  =  \ju\i//R'8ht  and  jTop  =  \r]\ysTop  allows  the 
relations  in  equations  (B.7)  through  (B.9)  to  be  changed  to: 


J 


Top 


W  -Slj 


■  Bottom 


■Right 


°Out 


So,,, 


■  Right  -Top  r,A  ci  A  ■ Left  -Bottom 

1  1  A  *-j  /  / 

-  +  Y~ - +Y  = - + - +  T-, - +  - 


\p\ex  W_ 


y 


a  a  juex 


(B.16) 


(B.17) 


(B.18) 


■'y 
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These  relations  can  then  be  solved  in  terms  of  the  outgoing  variables  as  shown  in 
chapter  four. 
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Appendix  C:  First  Spatial  Moment  Methods  Derivation  in  XY- 

Geometry 

Equation  Section  3 

This  appendix  contains  the  complete  derivation  for  the  first  spatial 
moment  method  presented  in  chapter  four  in  a  general  form.  Starting  with  the 
current  equations  using  sub-matrices  the  cell  face  equations  are: 


7 R  —  ;  L,  ,  ts-kk  ■ «  I  °  /3V 

JOut  ~  *-01  Jin  +  K0/  Jin  +  K0/  Jin  +  KOZ  Jin  +  *-O0UIn  +  *-O0Uln 


-RL  ~~:L 


-RRlR 


-RT  ~~:T 


RB~-B 


-  rl  ; 


+  KSd»  +  K«pf,  +K*mSj  +K*sa. 


057 


(C.l) 


~^L  t t-LL^L  .  -wt-LR-^R  .  vLT-T  .  1/ LB B  .  -w^LL  .  -wt-LR 

JOut  ~  KoiJln+KoiJin+KoiJin+KoiJin+KO0Vln  +  K 00Oln 


+KLJee)„  +K  L0‘e8-,„  +  K  LOSAs*  +k^s 


-LB  ] 


X  ,  \rL  o7  .  \rL  r'A 
+  iS-oSY^  +  *-OEAh  ’ 


(C.2) 


lout  =^Ollln  +KoiJ?n  +KoiJln  ^Ollln  +  Kg*?B 


+  Kga  +  +K^  +K  tosxSX  +K  tosySY  +K  toeaEa, 


(C.3) 


and 


(C.4) 


In 


-5  _  j^BL-'L  .  VBR~R  .  ^BT  ~~-T  .  ^BB~B  .  ^BL~AL  ,  VBR~AR  , 

JOut  ~  KOlJln  +  KOZ  Jin  +  KOZ  Jin  +  KOZ  Jin  +  K06?^"  +  KO0^»  + 

IS.  Off  a  In  +  i^oeVln  +  *-OSA^  +  *-OSX^  +  *-OSY^>  +  J^0£zfii  • 

addition  to  the  outgoing  currents,  the  outgoing  edge  distributions  are: 


fRL-.L  .  VRR~R  .  VRT~7 
U  Out  -  IS- 01  Jin  +  IV< 9/  Jin  +*<-01  J J„ 


RT  ~jT  ijsRB~jB.jrRLs\R 
-01  Jin  +  *-0I  Jin  +*-00  U J"  +*-00  & h 


+  K %0Tm  +  Kg*?B  +K  *saSa  +  KrsxSx  +  KrsySy  +  KreaE/ 

tout  =  K#7*  +  K#7*  +  Kg’. 71  +  KgJjJ  +  Kg  #B  +  Kg  *£ 

+  K##  Oin  +  +  K^-y^S  +  K 0EAE 


aT  _  i/TL-L  .  -it-TR  --R  .  VTT~T  .  VTB~B  ,  vTLnL  .  VTR~AK 
OOut  -  *<-01  Jin  +K0lJln  +K0lJln  +K0I  Jin  +*-000  In  +*-9edln 

-r-r-7  td-B  w^T  ~A  -x  ~Y  ,VT  LA 

L 0SY 0  +  *^0Ea1-  > 


+  k£S  +Kg^f„  +kSm5^  +kJ5X5x  +KjS7i 


(C.5) 


(C.6) 


(C.7) 
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and 


n  _  \{BL  ~jL  \  rs  BR  ~jR  i  t^BT  ~jT  .  t^BB  ~jB  .  j^BL  ~r\  _i_  \£  /i 

U Out  ~  IV 0/  Jln  +  JV,9/  JIn  +  IV 0/  JIn  +  JV#/  JIn  +f^go  v In  V In 

+  k| »Td»  +K^§f„  +K?mS-<  +k|sj,Sa'  +K|sl,Sy  +kIeaEa. 


The  cell  values  are: 


¥  -  K-x/j'In  +  KA7.//»  +  ^XlJln  +  ^XlJhi  +  K  V0 

+  +  +  KX&4S^  +  KxsxSA  +KX5FS}  +  KxeaEA 

¥  -  K Ylhn  +  ^YlJhi  +  K Yijjn  +  K  jlhn  +  Ky0Qln  +  K ygOln 

+  Y^yqOui  +  K.yq@!h  +  K-ysaS  +Kyyj^iS<  +KyyyiS  +  Ky^T  , 


S+kJ/75+k^+k^ 

+  K^X  +  K^5y  +  K^, 


+  KyffO^n  +  Kyfft 


(C.8) 


/  =  «47/«  +KiJ£  +Kr4/7;;  +  K*J*  +K^  +  K^ 

+  +K^  +K^  +  K^XSX  +  K^5r  +  KaeaEa, 

hi  In  +  ^XlJhi  +  K XI  jin  +  ^XlJhi  +  K X0 &I»  +  KX(9^« 

TXoO  in  +  KBx0Q[n  +  KxsaSa  +  KxsxSx  +  K  xsyS }  +^xea^A-> 


cA  v  ~A 
S  =£s¥  , 

S  =£s¥  , 

—Y  —Y 

and  S  =£sy/  . 

Equations  (C.12)  can  be  substituted  in  equation  (C.9)  to  eliminate 
average  scattering  source  which  gives. 


/  =(I-K  ASA  Y.S)~\K,]ln  +  42+45  +^Allf„  +-K^d»  +KRm(>1 

+  KtA00i>,  +  Kb40&i„  +  KasxSx  -  k  AS¥S ''  +K aeaEa], 

Letting  L4  =(I-K^S5)  1  equation  (C.15)  with  (C.13)  can  be  substituted  in 


(C.9) 


(C.10) 


(C.ll) 


(C.12) 


(C.13) 


(C.14) 


(C.15) 


equation  (C.10)  to  eliminate  the  x  moment  scattering  sources.  This  gives: 
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-X 


¥  -(I^i^xsx+KxsjLjKASX^s)  *[(KX/  +  ^XSA^sLA^AI)jji,+ 


(Kxi  +  K XSA  'Es  ^A^Al)jbi  +  (KA7  +  K XSA  ^ S  ^A^Al)jln  + 


~T 


(K xi  +  Kx&4  £ s  LAKB4j)j%  +  (Kx0  +  K xsa  £ 5  L^K^)(9/n  + 
(KX6>  +  KX&4  ’E'S  ^A^AdWln  +  (Kx6,  +  K  XSA  Zs  L^K^)^/,,  + 


(KV  ^  XS’/I  ^5  +  (Kxst  +  ^  XS’/I  ^  S’  L  A  ^  AS  y)S‘  + 


rB 


(C.16) 


(Xx£A  +  K XSA  I.V  LAKAEA)EA  ]. 


Letting  Lx  =  (I-(KXSi,  +  KAS4L  ^K^SY)ZS) _1  both  the  average  and  x 
moment  angular  flux,  equations  (C.15)  and  (C.16),  are  substituted  into  the  y 
moment  of  the  angular  flux,  equation  (C.ll),  with  equation  (C.14)  to  eliminate 
the  y  moment  scattering  source.  This  gives: 


¥  ~  (I  ~  (KKST  +  K  YSA  Z S  L.4  (K.4SF  +  K ASX  £ S  LX^XSY )  + 

K FSX  ^5  Lx(Kx5F  +  Kx&4  ^5  L^K^y))^) 

[(Kj7  +  Kf&4  Zs  L4  (Ki  +  K^5jy  Zs  LxKir )  +  kfsx  lx  (Ki 7  +  £,s  L.4K  i/  ))ii)  + 

(Ky/  +  Ky&4  Zs  (Kj,  +  K^jy  Zs  L  xk|j)  +  KF5X  Zs  Lx  (Kx/  +  Kx&4  Zs  L4K^7))y'^  + 
(Ky/  +  K  YSA  £  S  ^4  ^7/  +  Z5  LxK^y )  +  KF5X  Zs  Ljy  (Kx/  +  Kx&4  Zs  L^K^7))yJ;  + 
(Ky/  +KF&4Z5L4(K^7  +K^sxZ5LxKf7)  +  Ky5XZ5Lx(Kx/  +Kx&4Z5L4K^/))y7®  + 

(Ky0  +  K  YSA  Zs  La  (KAff  +  K ASx  Zs  LXKi"6»)  +  KF5X  Zs  Lx (KX6»  +  KXS A  Zs  LAKA0))0I„  + 


(K yo  +  KYsa  ~Zs  La (Ka0  +  KASX  Zs  LjyKxtf)  +  KF5X  Zs  Lx (KX(9  +  Kx&4  Zs  L + 
(kf<?  +Ky&4  Z5  L^(K^  +  ^ASX  ^5LXKX6»)  +  K  FSX  ^5Lx(KX6»  +  ^XSA  'E'S  L^TAQ))dln  + 


R 


(K Y0  +K YSA  ^ S  LA(Ka0  +  ^  ASX  £sLXKX6>)  +  K  YSX  ^sLx(^xe  +  L  XSA  'Z s  LaK Ag))01n  + 
O^YEA  +  ^  YSA  'ZjS  LA  O^AEA  +  ^  ASX  ^x^XEA  )  + 

KF5X  Lx(K^E4  +KX&4  ^5  L.4K.4£4))^j4]- 

(C.17) 


B 
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Let: 


Ly  -  (I  +  K ysa  'Es  La  (K asy  +  K asx  £ s  Lx^xsy) 

+  Kysx  S S  Lx  (K XSY  +  K  XSA  L aL- asy))^s) 

and  writing  the  following  for  the  first  term  in  equation  (C.17): 

mYI  =  Ly  (Ky/  +  Kf&4  Zs  La  (Kl4I  +  KASX  Zs  LxK lxi  )  + 

kfsx  'E'S  Lx(Kxj  +  Kxsa  Z5  L^7)). 

This  same  convention  is  followed  for  the  remaining  terms  and  allows  the  updated 
equation: 

n/ —  m vr  i t..  "I-  in irT  i /  -l-  m m  / r„  “I-  m w  /  r„  "I-  m nf) in  "l-  m  v /. f) jn 

(C.18) 

Equation  (C.18)  is  substituted  into  equation  (C.16)  to  eliminate  the  y  moment  of 


Yy  -  mYlJln  +  mYlJ In  +mYlJln  +mYlJln  +mY0& In  +  myffO/ii 

—*T  —*  B  —*■ 

+  m  Yffdln  +m  Yffdin  +m  yeaEA- 


the  angular  flux: 


XSY  +  L  XSA  £ S  La^ASY 


VX  =LA-[(Kij+KmZ.sL<K^+(K 

A  KA<^  X',  l-  Vi/  +  IK  ..,,  +  K  XS/I  La^ASY 
(Kj/+KmXsL4K^  +  (KXSY  +  K XSA  'Es  LAK ASY 


)Z s  mYl)Jlii  + 

)Z5my/)y'/f,  + 

)^SmYl)Jln  + 

^5  L.4KXSr)£,S  mYl)jin  + 


(Kaz  +  KXX4 


^s  L^A1  +  (KX57  +  kx&4 


(K^  +  K^ZsL4K^  +  (K 


XSY  +  K X&4  ^5 


)Z5  + 

(kJ,+  K  X&4  'Es  L^R4q  +  (KX5F  +  K XSA  'E'S  La^aSY^S  myoWln  + 

(  K  y  +  K  vv^  Zs  l-  Viy  +  '  ^  i;s  j  +  K  ,,.  ,  Ys  I.  in YffWln  + 


(Kjo  +  K 

XSA  £ s  La^a#  +  (K^sr  +  Kxsa  ^s  LaKasy)Hs  mYo)Oin  + 


B 


(K 


XEA  +  K  XSA  'E'S  LAK  AEA 


+  (K 


xsy  +KXSAT,sLAKASY)'LsmYEA)E‘  ]. 


(C.19) 


Let: 
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mA7  +K^ Zs +  (Kxsy  +  KXS4 Z5 mjj)  and  follow 

the  convention  is  for  the  remaining  terms.  This  allows  equation  (C.19)  to  be 
written: 


¥  -  m XI  jin  +  mX/  j/ii  +  ^XlJln  +  m XI Jin  +  m’x0 & +  mX0 

m  — >  J1  r)  — ►  A 

+  xQ@ In  "*"m XEA ^  • 

Both  equations  (C.18)  and  (C.20)  can  be  substituted  into  back  into  equation 


(C.20) 


(C.15)  to  eliminate  the  higher  order  scattering  moments.  The  average  angular 
flux  is: 


¥  ~  (L^K LAi  +  Kasx  Y,s  mA7  +  Z5  mYl)j/n  + 

(laka/  +  KAsx  £ s  mxi  +  KASY  Z 5  m  yi  )7/«  + 


(L.4KA7  +  KASx  Z s  mA7  +  K^5F  ^5  m  YI  )jjn  + 
(L^K^7  +  Kasx  Z 5  +  K^5F  Zs  m%j)jfn  + 

K.45X  £ 5  mX0  +  K^5F  mYe)d In  + 

(L^KA6» +  +  K^5F  Zs  niF6,)0/„  + 

(L^k^6»  +  k.4SV  +  Kasy  Z  5  mF6,)^7,;  + 

+  k.45X  ^5  vaBxe  +  K^5F  Z^  mf  e)din 


(L4K 


A^AEA  +  mXEA  +  KASY  S-,s 


£ SmYEA)E 


(C.21) 


Let  =  (L  jK^  +  Km,  Zv  m^7  +  K  /fSY  Zv  m |7 ) ,  and  follow  the  same 
convention  for  the  remaining  terms.  The  average  angular  flux  is: 


¥'  =  ™LAlJln  +  mAlJb:  +  mAlJln  +mAlJln  +  m70^«  +  m  Ad#In 


. R  AR 


T  lT 


B  AB 


. R 


* R 


+  +mBjn0i„  +m 


*5 


(C.22) 


M6> 


.40 c 


lAEA1 
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Now  equations  (C.22),  (C.20)  and  (C.18)  with  equations  (C.12)  through  (C.14) 
can  be  substituted  into  equation  (C.l)  to  eliminate  the  scattering  sources.  This 


gives: 


JOut  =  (K0/  +  K0&4  mU/  +  mA 7  +  KOSF  mYl)Jln  + 

(K of  +  Kosi  Us  +  Kqsx  mXI  +  mYl)Jbi  + 

(Kof  +  K OSA  £ 5  mU7  +  KOSX  £ 5  mX7  +  KOSF  mYl)jln  + 

Es  mU7  +  ^O.S'V  ^5  mf7  +  KpsF  mYl)jln  + 

(KOtf  +  K0&4  £ 5  mU6>  +  KOSX  ^5  mXi9  +  KOSF  ^5  mF<9 ln  + 
(KOtf  +  K0&4  ^5  mU<9  +  KOSX  mJF<?  +  KOSF  ^5  mYoWln  + 
(K06>  +  K0&4  ^5  mU<9  +  KOSX  mA7?  +  KOSF  ^5  mF^)^/«  + 
(KO0  +  K0&4  Ss  m^6»  +  KOSX  ^.S  mX<?  +  KOSF  mK6»)^/«  + 
(K0£4  +  KOS4  ^ .S'  m  .1 EA  +  KO.SV  ^s  m  XEA  +  KOSF  ^.S  m  YEA  )EA . 


(C.23) 


Let  ihqj  —  (K oi  ~t  ^^osa  ^s  ^  41  ^ os. v  Ss  nij^j  ^ctSF  Ss  and  use  the  same 


- RL 


convention  for  the  remaining  terms  in  equation  (C.23).  The  same  process  used  to 
produce  equation  (C.23)  is  applied  to  equations  (C.2)  through  (C.8).  The  result  is 
the  equations  for  the  sub-matrices  presented  in  chapter  four: 


lout  =  mOlJln  +  mOli‘!  +  mOI  jin  +  ™  Oil  In  +  mO0° +  m06°i> 
+  mO0  @In  +  mO0  @ln  +  mOEA E  ’ 


RLlL  ,  „ RRlR  ,  _ RT  ~~-T  ,  _ RB~;B  ,  __RL 


RR 


•R 


(C.24) 


~~-L  _W^LL-L  .  ^LR~R  .  ^LT~T  ,  ^LB~B  .  77  2»Z  .  ^LRnR 

JOut  ~  m  Of  Jin  +mOlJln  +mOlJln  +mOlJln  +  m()0& In  +  m0i9  6 In 

T  t  ~*T  r  n  -*  B  r  —■  a 

+  moe6In+mod6In+mOEAE'  , 


(C.25) 


7F  -mTL-;L  .TR-jR  .TT-jT  .TB-.B  .TLa1;  ,mTR  nR 

JOut  ~mOlJln+mOlJln  +  mOI J In  +mOlJln  +  mO0 “in  +  m()0t>l„ 
TT  TR  T  A 

+  mnn^ln  +  mnf)9ln  +  moEA E  ■> 


(C.26) 
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(C.27) 


-■B  _n%BL~L  .  ^BR~R  .  BT  ~-T  .  BB~B  ,  m BLnL  ,  m BR7>R 
J  Out  ~mOI  Jin  +  mOI  Jin  +  mOI  Jin  +  m()I  Jin  +  mOB@In  +  HI OO & In 

+  m 06  mOO  +  mOEA ^  ’ 


nR  _  mRL  ~~-L  .  ^RR~R  .  RT  ~T  ,  ^RB~B  .  mRL7>L  ,  m RRnR 
0 out  -  m  ei  j,n  +  m9I  jIn  +  m9I  jIn  +  m9I  j ,n  +  m99  6 in  +  mee  6 in 


+  mil  9}n  +  mil  6  In  +  m  nFAEA , 


_RB 


~B 


R 


/)  ___  1j I^j  •  Ij  *  __ _  JLI\.  •  A  .  __ _  jL/  J.  •  i  ■  __ _  J^D  *  D  ■  __ _  j /i  .  __ _  AA  /) 

0 out  -  m qi  j ,n  +m9I  jIn  +  m0,  j,n  +  m9I  jIn  +  mee  9 in  +  m99  6 in 

—*T  —*  B  —* 

+  moe  Oln  +  mee  Gin  +  mdEA^A  ’ 


_LR  --R 


LT  --T 


LB~B 


_LL\ 


.LR 


* R 


aT  -„JL~-L  ^™TR1R  ^™TT1T  ™trqR 

Oout  -m61jIn  +  m9Ijln  +  m0[  jIn  +  m9Ijln  +  moeOin  +meeOin 


and 


0 Out  -  me[  j,n  +  m0[  jIn  +mei  jIn  +mdI  jIu  +mee  9 in  +  meo  9 1„ 


.brar 


BTaT 


BBaB 


. bl ; 


. BR 


~R 


nr  •  / 

+  m^  9 in  + 


mBBdl  +  mBEAEA. 


(C.28) 


(C.29) 


(C.30) 


(C.31) 
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Appendix  D:  First  Spatial  Moment  Methods  Current  equations  in 

XY-Geometry 


Equation  Section  4 

Linear  Characteristic  Equations 

In  this  section,  the  relations  for  the  current  representation  of  the  linear 
characteristic  method  are  developed.  As  with  the  zeroth  spatial  moment 
methods,  the  usual  representation  found  is  in  terms  of  the  angular  flux.  In  the 
rectangular  cell  as  shown  in  figure  D.l,  the  equations  for  the  outgoing  angular 
fluxes  and  edge  spatial  moments  ij/,op  ,  x//nght ,  6top  and  0right  in  terms  of  the 
incoming  angular  fluxes  and  edge  spatial  moments  y/ho"om  .  y/h'n .  Qbottom  an(] 

9kft  scattering  within  the  cell  SA ,  Sx  and  S}  ,and  emissions,  EA  ,  are  previously 
derived  and  presented  in  the  literature.  (12:  23) 

Top 


Left 


Right 


Bottom 

Figure  D.l.  Rectangular  cell  for  first  spatial  moment  methods.  Cell 
shows  problem  variables  used  for  solving  the  discrete  ordinates  equations. 
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For  the  rectangular  cell  as  shown  in  figure  D.l,  the  equations  for  the 
outgoing  angular  fluxes  in  terms  of  the  incoming  angular  fluxes  and  edge 
distributions,  scattering  within  the  cell  and  emissions  are: 


y/Top  =  aM()(£y)y/Lefi  +(\-a)e  eyy/Bollom  + 

a[2Mx  (ey )  - M0 (£y )}0Left  -  a(\  -  a)e~£}  0Bottom  + 

fj[(l -a)M0(£y)  +  aMl(£y)]SA  + 

77  y  y 

^  «[-(!  -  a)M0  (£y)  +  (\-  2  a)Mx  (ey)  +  ocM2  (£y  )]S A  + 
||  [-(1  -  oc)M0  (£y)  +  (  2-  3  a)Mx  (£y )  +  aM2  (£y  )]S}  + 

77  y  y 

y/Right  =  M0  (£y  )YBottom  +  [(1  -  2  a)M0  (£y )  +  2  aMx  (£y  )pBoao/”  + 
^-M1(£y)SA  + 

77 

-  2  a)  M ,  (£)  +  2aM2(£)]Sx  + 

— [M^-M^S7  +^Ml(£y)EA, 

Tj  y  y  77 
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0Top  =  34(2ar  - 1  )M0 (4 )  -  2 aMx (£y)]ysLeft  +  3a(\  -  a)e  £y vBottom  + 

3a[(l  -  2  a)M0  (4  )  +  (6a-  2)MX  (4 )  -  4  ccM2  (ey  )]0Left  +  (l-3a  +  2a2  )e~£y  0Bottom  + 
— 34(1  -  oc)Mq  (e  , )  +  (2a- 1  )MX  (4  )  -  aM2  (4  )]£"*  + 

71 


^  [(1  -  3a  +  2a3  )M0  (e  v )  +  (3a  - 

71  y 


-  6a 3  )Mj  (4  )  +  6a3 M 2  (£y )  -  2 a3M3  (e  )]S x  + 


(D-3) 


Ay 

n 


34-(l  -  4A/0  (4 )  +  (3  -  (4 )  -  (2  -  5a)M2  (4 )  -  2«M3  (4  )]5,}  + 


— [(1  -  «)M0  (4 )  +  aM  1  (4)]^ , 

71 


0Right  =3[M0(£y)-2M1(£y)]i/7Bottom  + 

3[(1  -  2a)M0(£y)  +  (6a  -  2)M1(£y)-  4aM2(£y)\0Bottom  + 

—  3[Ml(£y)-M2(£y)]SA  + 
r)  >  y 


4 v 


M 


Av 


M 


Ay 

v 


3[(1  -  2a)M1  (4 )  -  (1  -  4a)M2  (4 )  -  2aM3  (4  )]5' A  + 
[-3MA4)  +  6M2(£y)-2M3(£y)]SY  + 
3[m1(4)-m2(4)]^, 


(D.4) 


y/A  -  aMx {£y)y/Lejl  +  [(1  -  a)M0 (4 )  +  aMx (4 )]^5oao"'  + 

43^  2{£y)-Mx  {£yj\0Left  +  a[-(l -a)M0(£y)  + (l  -  2a)M T  (4 )  +  aM  2  (4  )]<95oMom  + 

^[(1-4M1(4)  +  «M2(4)]^  + 

|/7| 

^  4~(1  -  i  (4 )  +  (1  -  2a)M2(£y )  +  aM3  (4  )]SA  +  (D-5) 

[-(1  -  «)Mi  (4 )  +  (1  -  2a)M2  (4 )  +  aM3  (4  )]S7  + 

—  [(1  -  «Wi  (4 )  +  ocM2  (£)~\Ea  , 

7/  y  y 
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y/x  =  3a[(2a  - 1  )MX  (ey )  -  2  aM2  (ey  )\y/Le^  + 

3a[(l  -  a)M0  (£y)  +  (2a  -  X)Mx  (ey)  -  aM2  (ey)\yBottom  + 

3a[(l  -  2 a)Mx  {ey)~  (4 a  - 1  )M2  (£y)~  2 aM 3  )]0ie-^  + 

[(1  -  3a  +  2a3)M0  (ey )  +  (3a  -  6a3  )MT  (^ )  +  6  a3M2  )  -  2  a3M3  )]6»SoWow  + 

^j3a[(l-a)M1(£-  )  +  (2a-l)M2(f )-aM3(e  )]5^  + 

77  •  (D.oj 

^  [( 1  -  3a  +  2a3  )MX  (ey )  +  (3a -  6a3  )M2  (ey )  +  6  a3M3  (^ )  -  2a3M4  (^  )]SX  + 

3a[-(l  -  a)Mx  (ey )  +  (2  -  3  a)M2  (£y )  -  (1  -  3  a)M3  (^ )  -  aM4  (fy  )]S}  + 

~3a[(l  -  a)Ml  (e  )  +  (2 a  -  1)M2  (f  )  -  aM3  (e  )]£' , 
rj  ■ 


=  3a[M !  (^ )  -  M  2(ey  )}y/Lel'  + 

3[(1  -  a)M0(f>,)  +  (3a  -  2)MX  (£y)~  2aM2  (£y)]i/fBottom  + 
a[-3M0  (ey )  +  6M,  (f}; )  -  2M2  (ey  )]0Left  + 

3a[-(l  -  a)M0  (^ )  +  (3  -  4a)MT  (ey )  +  (5a -  2  )M2  (ey )  -  2aM3  (ey  )]0Boltom  + 


Av 

V 


Ay 

v 


Ay 

n 


A y 


N 


3[(1-  a)Mj  (£’>,)  +  (2a-  1)M2  (^ )  -  aM3  (^  )]5y4  + 

3a[-(l  -  a)M1(f3,)  +  (2-3  a)M2(ey)~  (1-3  a)M3(ey)  -  aM4(sy)]Sx  + 
[-3(1  -  a)Mx  (ey )  +  (6  -  9  a)M2  (fy )  -  (2  -  8  a)M3  (^ )  -  2aM4  (^  )]5,}  + 
3[(  1  -  a)Mx  (£y)  +  (2a-  1)M2  (^ )  -  aM3  (^  )]F 4 . 


(D.7) 


For  these  relations,  /u  and  //  are  the  direction  cosines  along  the  x  and  y 


axis  respectively  from  the  angular  quadrature,^  = 


o-Ay 


M 


is  the  optical  thickness  in 
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the  y  direction,  a  =  —  is  a  parameter  for  the  cell,  and  M0(fv)  ,  M2(£y) , 

ex 

M3(£y)  ,  and  M4(£v)  are  the  exponential  moment  functions  (11:  27). 

The  definition  for  the  currents  jRlght  =  \ju\y/Rlght ,  and  jTop  =  \ri\y/Top  axvd 
similarly  for  the  left  and  bottom,  allows  the  transition  from  angular  fluxes  to  a 
current  representation.  For  the  first  spatial  moment  methods  however,  the  edge 
spatial  moments  must  also  be  transformed,  as  they  now  represent  the  spatial 
moment  of  the  current,  not  the  angular  flux.  This  is  done  in  the  same  manner  as 
the  currents  QRlght  =  | ju.\0Rlght  and  0Top  =\r/\dTop ,  although  the  notation  is  the  same 
for  both.  The  equations  (D.l)  through  (D.7)  in  a  current  representation  are: 

jTop  =!LaM  0(£. v)jLeft  +  (\-a)e~£y  jBottom  + 

M  ■  ' 

a[2Ml (£y ) - M0 (£y )]0Left  - a  (l  - a)e~£y 0Bottom  + 

A 

Ay[(\-a)MQ(£y)  +  aMl(£y)]SA  +  (D.8) 

Aya[-(  1  -  a)M0  (£y )  +  (1  -  2  a)Mx  (£y )  +  ocM2  (£y  )]5,A  + 

Ay[-(1  -  a)M0  (£y )  +  (2  -  3  a)M1  (£y )  +  aM2  (£y  )].S'y  + 

Av[(l  -  a)M0  (£y )  +  aM1  (£y  )]EA , 


■Right  =tLM  <e  \j Bottom  +lL[(\-2a)MQ{£v)  +  2aMl(£^}eBottom  + 
77  /7 

^^MX{£V)SA  + 
rj  7 

— [(l-2a)M1(£y)+2aM2(£y)]Sx  + 

rj 

^[M^-M^y^  +^MX(£y)EA, 

1)  ^  ^  1)  ^ 


(D.9) 
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eTop  =  —3a[(2a-l)M0(£y)-2aMx(£y)]jLeft  +  3a(l-a)e  ey jBottom  + 

A 

^-3a[(l  -  2  a)M0(sy)  +  {6a-2)Mx(ey)-4aM2{£y)\dLeft  + 

(]-3a+2a2)e~£>’6Bollom  + 

ky3a\{\-a)M()(£y)  +  (2cc-l)Mx(£y)-(xM2(£y)\SA  +  (D.10) 

Av[(  \-3a  +  2a3  )M0  (£y )  +  (3  a  -  6  a3  )MX  (£y)  +  6a3  M2  (£y )  -  2a3  M3  (ey  )]S A  + 
Ay3flf[-(1  -  oc)Mq  (£y )  +  (3  -  4a)Mx  (£y )  -  (2  -  5  a)M2  (£y )  -  2aM3  (£y  )]A}  + 

Av[(l  -  a)M0  (£y )  +  aMx  (£y)\EA , 

eRight  =  3 [M0  (£  )  -  2 Mx  (£v  )]jBottom  + 

V 

^3[(\-2a)M0(£y)  +  (6a-2)Mx(£y)-4aM2(£y)WBottom  + 

—3 [Mx(£)-M2(£)\Sa  + 
n 

A  (D'U) 

~^"3[(1  -  2  a)Mx(£y)-  (1  -  4  a)M2(£y)-  2  aM3(£y)]Sx  + 

^y[-3Mx  (£y )  +  6M2(£y)~  2  M3(£y)]SY  + 

— 3[Mx(£y)-M2(£y)]EA , 
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¥A  =  -aMx(ey)jLeft  +-\_{\-a)MQ{ey)  +  aMx(£y)]jBottom  + 

//  -  rj 

-(AM2{ey)-Mx(eyWLeft  + 

M 

|[-(1  -  a)M0(ey )  +  (1  -  2a)Mx{sy )  +  aM2  (£y)]0Bollom  + 

— [(1  -a)M1(ey)  +  aM2(ey)]SA  +  (D.12) 

—[-(1  -  a)M  i(£y)  +  (1  -  2a)M2(£)  +  aM3(£y)]Sx  + 

rj  - 

— [-(1  -  oc)Mx  (£  )  +  (1  -  2 a)M2  (£  )  +  aM3  (£  )]S}  + 

rj  -  y 

— [(1  -CC)MX(£  )  +  aM2(£  )]Ea, 

rj  >  > 

¥x  =-3a[(2a-\)Mx(£y)-2aM2(£y)]jLefi  + 

// 

ha[{\-a)M0{£y)  +  {2a-\)Mx{£y)-aM2{£y)]jBottom^ 

Tj 

-  3a[(l  -  2  a)Mx  (£v)~  (4a  - 1  )M2  (£v)~  2  aM3  (£v  )]0Left  + 

//  •  y  $ 

—[(1  —  3 a  +  2 a3 )M0  (£  )  +  (3a  -  6a 3 )MX  (£  )  + 

rj  y  y 

6aiM1  (£y )  -  2o?M3  (£y  )]0Bottom  + 

— [(1  -  a)Mx(£y ) +  (2a  - 1  )M2(e)-  aM3(e)]SA  + 

n 

— [(1  -  3a + 2  a3  )MX  (£y )  +  (3  a-  6a 3  )M2  (£y )  +  6  a3M3  (£y )  -  2  «3M4  )],S  A  + 

^[-(1  -  )  +  (2  -  3  a)M2  (£y)  -  (1  -  3  a)M3(£y)  -  aM4  (£y)]SY  + 

—[(\-a)Mx(£y)  +  (2a-l)M2(£y)-aM3(£y)]EA,  (D-13) 
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VY  =  —\Ml(£y)-M2(£y)\jLeft  + 

M 

-[(1  -  a)M0(e)  +  (3a-2)Ml  (e)  -  2aM2(e)\jBoUom  + 
i) 

a[-3M0 (ey )  +  6Mj (£y)~  2 M2  (£y  )]0Left  + 

^[-(1  -  a)M{)(£y)  +  (3  -  4a)M,  (£y)  +  (5a-2)M2(£y)  -  2aM3(£y)]0Bottom  + 
^[(\-a)Ml(£y)  +  (2a-\)M2(£y)-aM3(£y)]SA  + 

^[-(1  -  a)Mx  (£y)  +  (2  -  3a)M2(£y)  -  (1  -  3a)M3(£y)  -  aM4(£y)\Sx  + 
—[-3(1  -  a)Ml  (£y)  +  (6-  9  a)M2  (£y)  -  (2  -  8  a)M3  (£y )  -  2  aM4  (£,>,)]5'1  + 
— [(1  -  a)Ml  (£y)  +  (2  a- 1  )M2  (£y )  -  aM3  (£y)~\EA . 


These  are  the  relations  used  in  chapter  four  for  substitution  into  the  first 
spatial  moment  method. 

Linear  Discontinuous  Equations 

In  this  section,  the  relations  for  the  current  representation  of  the  linear 
discontinuous  method  are  developed.  As  with  the  zeroth  spatial  moment 
methods,  the  usual  representation  found  is  in  terms  of  the  angular  flux.  In  the 
rectangular  cell  as  shown  in  figure  D.l,  the  equations  for  the  outgoing  angular 
fluxes  and  edge  spatial  moments  y/top  ,  y/"s'u  ,  6top  and  0""1"  in  terms  of  the 
incoming  angular  fluxes  and  edge  spatial  moments  i//ho"om  .  y/left .  Qbo,,om  an(] 

0,eft  scattering  within  the  cell  SA ,  Sx  ,  SY  ,and  emissions  EA  ,  are  previously 
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derived  and  presented  in  the  literature  (4:  289-290).  The  linear  discontinuous 
relations  for  the  edge  values  are: 


l//Right  =y/A+lf/x, 

(D.15) 

y/Top  =  y/A  +  i//y , 

(D.16) 

eRight  =  y/Y, 

(D.17) 

and 

dTop=y/x. 

(D.18) 

These  can  be  substituted  into  the  zeroth,  x  and  y  moment  cell  balance 
equations: 


a{y/Right  -  VLe,t)  +  (y/Top  -  y/BoUom)  +  £wA  =sA^, 

V 

(D.19) 

3 a{y/Right  +  VLeft  -lYi)  +  (0T°P  ~ 0Bottom)  +  £vy/X  =  SX^, 

y  7) 

(D.20) 

3 {x/fTop  +  VBott°m  -  )  +  cc{dRight  -  0Left)  +  £xm/y  =  SY^, 

y  ij 

(D.21) 

where  the  following  relations  are  defined: 

,  3  3er 

Cl  —  1  +  OC  +  £.,  H - 1 - , 

4  +  £y  1  +  3  OC  +  £y 

b  =  4  +  £y,  (D.22) 

c  =  1  +  3  a+£y. 

After  some  algebra,  the  cell  values  are: 


r 


AySA  +  y/Bo“om  +  awLeft  + 
V 


3  y/ 


Bottom 


_  0Left  _fysY  a{3y/Left  _  & Bottom  _^X 


V 


-+- 


V 


(D.23) 


184 


Ay  O  X  ,  a Bottom 


¥ 


n 


SA  +  gUottom  +  2a(i^A  _  yLeft^ 


(D.24) 


and 


¥Y  = 


Ay 

v 


sY +eLeft +\ya -¥Bottom) 


(D.25) 


The  desired  relations  in  the  angular  flux  relation  are  found  using  equations 
(D.15)  through  (D.18)  with  equations  (D.23)  through  (D.25): 


y Top  _  |-(3  +  b)a  +  3(3  +  b)a^Left  +  (9  +  6b-3ab  +  b~)  ^ Bottom  +  (-3 +  (<?-!)£)  ^Left 


ab 


abc 


ab1 


ab 1 


a(b  +  3)  e Bottom  +  Ay(h  +  3)  _  Aya(b  +  3)  ^  +  Ay(-3  +  (a-\)b)  ^y  + 


abc 

Ay(h  +  3)  A 
ab  rj 


ib\r/\ 


ahc|^| 


ab2  rj 


(D.26) 


yjRight  _  r  3a(c  ac  +  3 a)  ^  a(c  +  3a)  -.  Left  _|_  (3  +  b)(c  +  3a) 

'  L  9  -I'  /  ' 

abc 


Bottom 


ac 


ac 


(c  +  3a) 
abc 


gLeft 


+ 


(ac-a(c  +  3a))  |  Ay(c  +  3a)^4  |  Ay(ac-a(c  +  3a))^x  _ 


ac 


ac\tj\ 


2  „ 
ac  rj 


(D.27) 


Av(c  +  3a)//  +  Ay{c  +  3a)  eA 


ahc|//| 


ac\rj\ 


qTop  _  j- 3a  _  3a(ac  3a)]^jjg^  +  3(3  +  h)a  ^ Bottom  _  3or  gLeft 


ac 


ac 


abc 


abc 


+ 


(ac-3a“) 


Ay3a  sa  |  Av(ac-3q-)^x  Ay3a  (  Ay3a 


ac\rj\  ac~  \r/\  abc\ii\ 

nRisht  r3ct  9a  9  +  3B  —  3ab 

9  8  =  [ —  + - w  + - 7 - ¥ 


ac\7l\ 


ab ' 


Bottom  .  (  3  +  ab) 
ab2 


gLefi  - 


3  a-e  Bottom  _j_  3 Ay  sa  3A ya  (  Ay(ab-3)  sy  {  3Ay 


ahc 


a 


%l 


ahc  r] 


ab2  r/ 


ab  r] 


ac~ 


0 


Bottom 


+ 


(D.28) 


(D.29) 
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¥ 


A  _  j- &  _j_  j yfLefi  _|_  (3  +  b)  Bottom _ ]_ gLeft  _  q 

a  ac  ab  ab  ac 


i  Bottom 


+ 


Ay  cA  A ya  ?X  a y  SY  |  Ay  ea 


hi 


ac  77 


ri \ij\ 


am 


x  =  _  3a(ac-3a)  Left  +  3(3  +  b)a 

1  \-  ?  7  • 

abc 


Bottom 


ac 


ac 


3  a 
abc 


(ac-3er) 


gLeft  \wv 


j Bottom 


ac' 


Ay3a  sa  |  Ay(ac-3or)  sx  Ay3a  |  &y3a  £A 


ac\m 


ac2  77 


a6c|//| 


act; 


¥ 


Y  _  j-30f  -^Left  _j_  (9  +  36  3a6)  ^ Bottom  _j_  (  3  +  ab)  gLeft 


"  ab  abc 


ab1 


ab1 


3  a 
abc 


e 


i  Bottom 


+ 


Av3  ^  Ay3or  ^  |  Ay(afr-3)  r  |  Ay3or^ 


(D.30) 


(D.31) 


(D.32) 


a6|^|  a6c|^|  a6“|^|  ac\f]\ 

The  definition  for  the  currents  jR,ght  =1  ju\y/Rlght ,  and  jTop  =  Yi\¥T°P  and  similarly 


for  the  left  and  bottom,  allows  the  transition  from  an  angular  flux  to  a  current 


representation.  As  with  the  LC  method,  the  edge  spatial  moments  must  also  be 


transformed,  as  they  now  represent  the  spatial  moment  of  the  current,  not  the 
angular  flux.  This  is  done  in  the  same  manner  as  the  currents  QRlghl  =  \fl\ QRlght , 
and  eTop  =  \t]\  0Top  although  the  notation  is  the  same  for  both.  The  equations  for 
the  outgoing  quantities  are: 


(D.33) 


abc 


ab 


abc 


■Right 


a(c2  +9a  +  c( 3  -  3a  +  3 a))  .Left  (3  +  b)(c  +  3 a)  |//|  ,Bottom  (c  +  3 a)  Left 
1  J  +  ,  LJ  J  ,  “  ' 


ac 


(ac  -  a{c  +  3a))  |//|  ^Bottom  +  Ay(c  +  3a)\/u\  ^  ^  Ay(ac  -  a(c  +  3a)) /u  c x 


ac 2  r/ 


abc\rj\ 

SA+- 


ac\rj\ 


ac2r/ 


abc 
SA  - 


(D.34) 


Av(c  +  3cf) |/y |  Y  |  Ay{c  +  3a)\fi\  A 


abc  7] 


ac\rj\ 


pTop  _  3a(  aC  +  (3  +  c)a)  //  .Ze//  3(3  +  .Bottom 

ac2 1//|  J  abc  J 

gLeft  _j_  ~  36y~  )  Bottom 

abc\ju\  ac 2 

Ay3a  s A  |  Ay(ac-3ar)  Ay3ar  y  |  Ay3a  rA 

ac  ac 2  abc  ac 


(D.35) 


_  3a(3  +  c)  ; iefi  _  3(  3  +  (a  1)/?)  //  ;Bottom  ^  (  3  +  ab)^Lef, 


QKignt  _ 


abc 


ab 2  r/ 


-j  +- 


ab2 


3a\ju\ 

abc\rj\ 


e 


Bottom  _|_  3Ay\ju\sA  3Aya\ju\sX^ 

ab  rj  abc  /] 


Ay(ab-3)\fl\  Y  3Ay|//|  4 

O  H - — — —±L  , 


a62|/;| 


7^7 |/7| 


_  C({3  +  C )  ,z.e//  (3  +  6)  ■  Bottom  1  pLeft 

V/  = - 1 — i-/  + - 1 — r  - 1 — r"  —  ■ 

ac|//|  a6|^|  a6|//| 

A y_sA  Aya  Ay  c;}-  |  Ay  ^ 


^  pBottom 


ac  7 7 


ac  77 


?6|//| 


a  77 


+ 


,,,x  3or(  ac  +  (3  +  c)or)  3(3  +  b)a  .Bottom 

V  = - 2^ - J  + - 


ac"  // 


abc  |/7| 


3or 


a6c  ju 


gLeft 


+ 


(ac-3a  ) 
ac2  ti 


6 


Bottom 


+ 


Ay3a  s A  |  Ay(ac-3a~)  sx  Av3a  y  |  A y3a  £A 


ac  77 


ac2  7] 


abc  |/7| 


ac  77 


(D.36) 


(D.37) 


(D.38) 
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Y  _  3or(3  +  c)  jLeft  |  (9  +  36  3ab)  .Bottom  ^  (  3  +  Clb )  gLeft  3 O'  ^ Bottom 

abc\ju\  ab2  \tj\  ab2\ju\  abc\fj\ 

Av3  ^  Ay3or  a-  ,  Ay(a6-3)  |  Ay3or  ^ 

aZ?|^|  aZ?c|^|  a/j2|^|  <2c|//| 

These  are  the  relationships  used  in  chapter  four. 


(D.39) 
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Appendix  E:  First  Spatial  Moment  Methods  Derivation  in  Slab 

Geometry 

Equation  Section  5 

This  appendix  contains  the  complete  derivation  for  the  first  spatial 
moment  method  presented  in  chapter  three  in  a  general  form.  The  first  spatial 
moment  methods  need  several  additional  equations  to  account  for  the 
contribution  to  the  flux  from  the  first  moment  of  the  scattering  source.  The 
relations  described  for  the  zeroth  spatial  moment  methods  in  chapter  three 


become: 

=K0/(?"  +IW-4  +Kosj,Sa'  +KoeaEa  ,  (E.l) 

VA  =  KAIp" +KasaSa  +KasxSx +KaeaEa  ,  (E.2) 

¥X  =^xi¥m +^-xsaSA +K-xsxSX +K-xeAE  ,  (E.3) 

SA=ZSWA,  (E.4) 

and 

Sx=Y.syx.  (E-5) 


Again.  K0/ ,  K()V( ,  K0&) ,  KAl ,  K  ASA ,  K  ,.i  V  ,  KAEA ,  K  A/ ,  K  .  .  . ,  K  Avv  ,  and 
K xea  represent  diagonal  matrices  of  transport  coefficients  that  define  the 
relations  of  the  inputs  of  a  cell  to  the  calculated  quantity.  For  example  KX1 
represents  the  contribution  to  the  first  moment  flux  from  the  incoming  flux  and 
is  the  scattering  matrix  described  in  chapter  three.  These  matrices  are  the 
sub-matrices  used  in  the  general  derivation  shown  in  chapter  three. 
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Equations  (E.4)  can  be  substituted  into  equation  (E.2)  to  eliminate  the 
average  scatter: 

r'  =(I-K^ZS)_1[K^“  +K ASXSX  +K aeaEa]  .  (E.6) 

As  a  shorthand  notation,  let  LA  =  (I-K ASA  Z^)-1  •  This  result  can  be  substituted 
into  equation  (E.3)  with  (E.5)  to  eliminate  both  the  x  moment  of  scatter  and  the 


average  scatter: 


Let: 


and 


¥X  -  ( 1  -  ( K  XSX  ^  S'  +  K  ASX  £  S'  L  ,4  K  ASX  ))  lx 

— -in 

[(Kx/  +  K ASX  Z5  L  4K AI)y/  + 

—*A 

(K-XEA  +  K ASX  £ S  LAK-aEa)E  ]• 


-  (I  -  (KXS X  'E'S  +^ASX  'E'S  Ia*asx  Zs))  1  > 

mXI  =  (KX/  +  KXSX  'Es  ) 


(E.7) 


(E.8) 

(E.9) 


mXEA  -  (KxEA  +  KasX  £ s  LaK-aea)  , 

are  used  for  equation  (E.7).  Equation  (E.7)  is  now  substituted  back  into 


(E.10) 


equation  (E.6)  to  eliminate  the  x  moment  of  scatter: 

\j/A  =  L4[(K4/  +K  asx  Us  ^XmXlW"  + 
(KAEA  +KASX 


Let: 


m  4I  ~  (KAI  +  KASX  'E'S  LXmX/)  > 


and 


(E.ll) 


(E.12) 
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mAEA  ~  (K-AEA  +  KASX  £ S  ^XmXEA)  > 


(E.13) 


are  used  for  equation  (E.ll).  Both  equations  (E.7)  and  (E.ll)  can  be  substituted 
in  equation  (E.l)  to  eliminate  the  scattering  terms.  This  gives  the  outgoing 
detailed  flow  for  a  cell  in  terms  of  the  incoming  detailed  flow  and  emissions  in  a 
cell: 


¥°Ul  ~  (K0/  +  K osa  £ s  L AmAi  +  K osx  'Es 
+(KOEa  +  K-OSA  'Es  ^AmAEA  +  K-OSX  L XmXEA  w. 


(E.14) 
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